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Abstract 

The magnetorheological (MR) damper behavior is 
described by a set of nonlinear differential 
equations, which are unsuitable in control 
applications. This paper illustrates the effectiveness 
of adaptive neuro-fuzzy inference system 
(ANFIS) in building both forward and inverse MR 
damper models. The forward model is obtained by 
training noisy data generated from simulation 
of the nonlinear differential equations with low pass 
filter. Input selection method is applied to find the 
most appropriate inputs to train the inverse 
model. Finally, the numerical simulation proves 
that the two models established based on ANFIS 
can track the MR damper force for different 
amplitude and frequency. In addition, we can 
predict the voltage required to be applied to 
produce a known force determined from control law 
with high performance. 

Keywords: ANFIS, Forward fuzzy model, Inverse fuzzy 
model, Magnetorheological damper. 

1. Introduction 

Magnetorheological damper is a semi-active device 
that requires minimal power (20-50 watt); this 
device can operate with battery eliminating the 
need for a large power source. Because the device 
forces are adjusted by varying the strength of the 
magnetic field, no mechanical valves are required 
making a highly reliable device [1]. 
MR dampers are used in various applications, which 
have been and continue to be developed. These 
dampers mainly used in industry such as heavy 
motor damping, vehicle suspension [2]-[7], seismic 
protection of building during earthquakes [1], [8]- 



[18], and the operator seat damping in the 
construction vehicle [19]. Other fields use MR 
damper is the military and commercial helicopter 
cockpit seat. Another application humanitarian 
is the prosthetic leg. Its decrease the shock 
delivered to the patient’s leg when he is jumping 
[ 20 ], [ 21 ]. 

Parametric models such as Bingham's model, 
extended Bingham, and Bouc-Wen, which 
numerically tractable and has been used extensively 
for modeling hysteretic systems [22]. Regardless of 
the advantage of 

this model, the force-velocity relation within the low 
velocity range does not closely reproduce the 
measurement results obtained from a typical MR 
damper [23]. So, based on extensive experiments, 
Spencer J. et al. [24], [25] proposed a modified 
Bouc-Wen model, which avoided most of the 
problems observed in the other models. This model 
could keep up all advantages exhibited by the 
Bouc-Wen model. The parameters associated with 
this model were found through matching the 
simulated data with experimental data obtained in a 
variety test. 

Several techniques are used for emulating non 
parametric models such as neural networks, fuzzy 
logic, neuro-fuzzy and polynomial models. Wang 
and Liao [26] proposed a forward MR damper 
model, which consists of a recurrent neural-network 
in which the output is delayed and fed back to the 
input layer. The Levenberg- Marquardt algorithm is 
the selected to be training method and the training 
data was obtained from the phenomenological 
model developed by Spencer J. et al. [24], [27]. The 
results show that the predicted damping 
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force approximates the target data except at low- 
frequency. However, as mentioned by the authors, 
these are still a preliminary work and further 
research is required before practical applications. 
[28]- [30] presented a forward neuro fuzzy 

model for large-scale and 30 ton MR damper. 

The reminder of this paper is organized as follows: 
Section (2) focuses on the model analysis of MR 
damper. Sections (3) describes and discusses 
the forward fuzzy model. While, section (4) 
describes and discusses inverse fuzzy 
model. Finally, the conclusion is presented in 
section (5). 



2. Model analysis of MR damper 



MR damper is consisting of the main 
hydraulic cylinder which, houses the piston, the 
magnetic circuit, and accumulator. The cylinder is 
filled wirh MR fluid which, composed of micron- 
sized magnetically polarizable particles dispersed in 
a carrier medium such as water, mineral or synthetic 
oil. Typically, it contains 20 to 40% by volume of 
relatively pure carbonyl iron with 3 to 5 microns in 
diameter [31]. MR fluid is normally a free flowing 
viscous fluid, but the presence of a magnetic field 
causes the particles to form chains and increase the 
fluid viscosity, until it becomes a semi-solid (Figure 
l.a). 




(a) (b) (c ) 



Figure l.a: Operation of MR fluid (a) no magnetic field (b) and (c) 
with increasing magnetic field. 



The modified Bouc-Wen model is shown in 
Figure. l.b. It is consisted of mechanical elements 
such as springs and dashpots to emulate the device 
behavior. 



y x 




Figure l.b: Modified Bouc-Wen Model. 



The model is described by the following seven 
nonlinear differential equations [25]: 



F =c l y + k l (x-x 0 ) (1) 

where 

F : The damping force. 

The viscous damping for force roll-off at low 
velocities. 

iq:The accumulator stiffness. 
x : The relative displacement. 

: The initial deflection of the accumulator gas 
spring. 

The internal displacement of the damper y and the 
evolutionary variable z are governed by the 
following equations: 

y= — - — {az +c o x +k o (x -y)} (2) 

+c i 

Reference [32] described the restoring forces with 
hysteresis in the following form: 

z =-y\x ~y \z " _1 \z \- P(x -y)z n +A(x -y) 

n even (3-a) 

z =-/\x -y\z n - P(x - y)|z"|+A(i - y) 

n odd (3-b) 

where 

C o : The viscous damping at large velocities. 

k o : The viscous damping for force roll-off at low 

velocities. 

X : The piston velocity. 

a : The Bouc-Wen parameter describing the MR 
fluid yield stress. 
n : Constant. 



The adjustment of hysteresis parameters a, p, and A 
determines the linearity in the unloading region as 
well as the transition smoothness from the pre-yield 
to the post-yield region. The parameters y, ft, A, n 
and k 1 are considered fixed and the parameters c 0 , Ci 
and a are assumed to be a function of the applied 
voltage u . 

a = a a + a b u (4) 

= C la + C lb U (5) 

= C oa +C ob U (6) 

If U c is the command voltage sent to the current 

driver, the dynamics involved in the MR fluid 
reaching rheological equilibrium is modeled by the 
first order filter. 

u=-r/(u-u c ) (7) 

where 7j is the time constant. 

Dyke [17] shows the response of the model is 
dominated by the input velocity and is extremely 
sensitive to change in velocity. This property causes 
the model behave erratically when there is a sudden 
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change in the velocity. Then a Low Pass Filter (LPF) is 
added to improve the ability of the model to deal 
with noisy displacement data. 

The dynamics of the MR damper with LPF becomes: 

y= — - — {az +c a v +k 0 (x -y)} (8) 

C o +C 1 

z =~y\> -y'\z n ~ l \z\- Piy -y')z n +(v -y') (9) 



w =A f w +B f x 
v =C f w 



( 10 ) 



A ~ 



0 



-ox 






co. 



B,= 



C, -[* 0] 



Where Ay , Z?y and Cy are system matrices of LPF, 
W is the state of the filter and v is the filtered 



velocity, CO c is the cut off frequency. The 14 

parameters which characteristic the MR damper are 
given in Table I. 



The ANFIS has five layers [33], in which node 
functions of the same layer have the same function 
type as described below: (Note that Oy denotes the 
output of the i th node in the j th layer). 

Layer 1: Every node i in this layer is an adaptive 
node with a note output defined by 

O u = (x ),i =1,2 or 

( 11 ) 

°u = B Bi _ 2 (y ). 1 = 3 > 4 

where x (or y) is the input to the node and A t (or B^) 
is the fuzzy set associated with this node 



/**,(*) = ■ 



i+ 



f \ 2 
x -c, 



V a i J 



( 12 ) 



Where, | a { ,b i ,c. } is the parameter set premise 
parameters. 



Table I: Parameters for MR-Damper 



Parameter 


Value 


Parameter 


Value 


COa 


21.0 N. sec/cm 


a a 


140 N/cm 


Cob 


3.50 N. sec/cm. V 


a b 


695 .sec/cm.V 


k 0 


46.9 N/cm 


P 


363 cm' 2 


Cla 


283 N.sec/cm 


P 


363 cm' 2 


Clb 


2.95 N.sec/cm.V 


A 


301 


ki 


5.00 N/cm 


n 


2 


X 0 


14.3 cm 


T] 


190 sec' 1 



3. Forward fuzzy model 

Neuro Fuzzy model is a non parametric model for 
emulating MR dampens behavior. It is faster than 
the mathematical model for keeping the error 
relatively small. The neuro fuzzy is an application of 
ANFIS. ANFIS model using neuro adaptive learning 
techniques, which are similar to those of neural 
networks, were originally presented by Jang [33], 
[34]. 



3.1Anfis Architecture 

The architecture of a two-input two-rule ANFIS, as 
example is shown in Figure 2. 

k y 




Figure 2: ANFIS architecture. 



*Layer 2: Every node in this layer is a fixed node 
labeled II, which multiplies the incoming 
signals and outputs the T-norm operator 
result, e.g. 

°i,i =w i = Ba, ( x ) x fi Bt (y ), i = 1, 2. (13) 

Each node output represents the firing strength of a 
rule. 

® Layer 3: Every node in this layer is a fixed node 
labeled N 



0 3,i = W , 



W j +W 2 



J 



U. 



(14) 

Outputs are called normalized firing strengths. 



■Layer 4: Every node i in this layer is an adaptive 
node with a node function 

°4,i =w if i =™i ( p t x +^iy +0 (i5) 

where Wj is the output of layer 3 and {/? ( ,4, ^} is 
the parameter set consequent parameters. 



1 Layer 5: The single node in this layer is labeled 27, 
which computes the overall output as the 
summation of incoming signals. 



Yx if 

°5A =YXifi = V 

Y w i 



( 16 ) 
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3.2 Training Forward model 

Given input/output data sets, ANFTS constructs 
Fuzzy Inference System (FIS) whose membership 
function parameters are adjusted using a back 
propagation algorithm. The size of input-output 
data must be large enough and cover all frequency 
and amplitude ranges to fine tuning the 
membership function. Three inputs corresponding 
parameters namely; displacement, velocity, voltage 
and the MR damper force was the output, which 
used for the ANF1S training and checking data to 
build the forward fuzzy model. The training and 
checking data are obtained from simulating the 
mathematical Equations (1), (4)-(10) into MATLAB 
program. The displacement inputs generated using 
chirp signals with amplitude 4 cm and frequency 
between 0-3 Hz. The voltage inputs generated 
using Gaussian white noise ranges from 0-3 V with 
frequency 0-3 Hz. Our simulation was based on off- 
line learning only. The data were collected for 20 
second and sampled at 1000 Hz, which generate 
20000 points of data. 

The odd sequence numbers were chosen to be the 
training data (10000 points) while the even 
sequence numbers were used as checking data 
(other 10000 points). The power spectrum density 
for input and output data are shown in Figure 3. 

500 
0 

-500 

_ 0 50 100 150 200 250 300 350 400 450 500 

'n' 

X 

S 100 

1 
3 



t -100 

I 0 50 100 150 200 250 300 350 400 450 500 

c. 

100 
0 

-100 

0 50 100 150 200 250 300 350 400 450 500 

Frequency(Hz) 

Figure 3: Power spectrum density of displacement, voltage and 
force 



(c)Force power spectrum density 




■ 1 i i i i L 



(b)Voltage power spectrum density 




(a)Displacement power spectrum density 

t 1 1 1 1 1 1 1 r 

* 

J I I I I I I I L 



3.3 Validation of forward model 

To validate the accuracy of the forward model, five sets 
of data are generated from mathematical equations to 
determine the target forces. The first and second sets are 
sinusoidal excitation with 3cm amplitude, two different 
frequency 2.5 Hz for data set I and 1 Hz for data set II 
and held for three constant voltage 0V, IV, and 2.25 V. 
Figures (4) and (5) show the comparison between time 
history of target and predicted forces and force-velocity 
relationship. It can be seen that the predicted force can 
track the target force very well for both 0V and 2.25 V 



but in the case of constant 1 V the hysteresis behavior 
predicted not coincide with the target. 

The other three validation data sets are chirping 
displacement with step, pulse, and sinusoidal voltage. 
For this validation sets the voltage time history is 
plotted in Figure (6). Figures (7)-(9) show the 
comparisons between the target and predicted forces 
versus time and velocity. The predicted forces and 
hysteresis behavior are matched excellently with the 
target for a varying voltage with time( the case of semi- 
active). The accuracy of the neuro fuzzy forward 
model is proven theoretically by using the error 
analysis which described by Spencer et al. [25]. The 
error analysis for all data sets are shown in Table n. 
E t ,E v , and E a are the normalized model error 

as a function of time, displacement and velocity. 
From this table, the results show that the force 
predicted by forward fuzzy model can predict the 
target force for all cases very well. 




Amplitude=3cm, Frequancy=lHz 




(b) 



Figure 4: Comparison between fuzzy and mathematical model for data 
set I: (a) force-time history ; (b) force-velocity relationship. 



Table II: Error norm for validation of the forward fuzzy model 



Data 

Sets 


Disp. 


Volt 


E, 


E v 


E a 


I 


3sin(57it) 


0 


0.0004 


0.0027 


0.0007 


I 


3sin(57it) 


2.25 


0.0002 


0.0017 


0.0004 


n 


3sin(27it) 


1 


0.00003 


0.0003 


0.0003 


n 


3sin(27it) 


2.25 


0.0003 


0.0012 


0.0002 


in 


Chirp(0-3)Hz 


step 


0.0002 


0.0021 


0.0005 


IV 


Chirp(0-3)Hz 


pulse 


0.0077 


0.0673 


0.008 


VI 


Chirp(0-3)Hz 


l+0.5sin7i;t 


0.0007 


0.0004 


0.0008 
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Table III: comparison between the error norm for the present and 
previous work [29] 



Disp. 


Volt 


E t 




A 




E a 




2sin47it 


0 


0.078 


0.037! 


0.271 


■ 


1.244 


§533 


2sin47it 


2.25 


0.027 


ImI 


0.098 


0.086 


0.541 


0.875 


0.5sin7i;t 


0 


0.013 


0.471 


0 m 4 


0.441 


0.001 


0.904 


0.5sin7i;t 


2.25 


0.011 


|.235 


0.014 


1.264 


0.001 


0.365 


GWN 

0-2Hz 


GWN 

0-2Hz 


o.ooi imH 


o 

o 


ffl 


o.ooi H 



0 The present work Q The previous work[29] 



A comparison between the present ANFIS forward 
model and forward model created by other [29] are 
shown in Table ffl. It is show that the values of 
calculated error norm are very small compared with 
the same set of data and the same time span. 



Step voltage 

3 1 1 1 1 1 1 1 1 

2 - 
1 ■ 

0 i i i i i i i i 

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 





Figure 6: Voltage time history for data sets III, IV, and VI. 



4. Inverse fuzzy model 

Control of the MR damper is very important for all 
applications. Because the MR damper force cannot 
be directly controlled, then we can control the 
applied voltage which, can be directly controlled to 
operate the MR damper. To achieve that, we must 
build an inverse model. It predicts the voltage 
required to be applied to the current driver to 
produce a desired force with known displacement 
and velocity. 



Amplitude=3cm, Frequancy=2.5Hz 





Amplilude-3cm, Frcqiiancy-2.5H£ 



-50 



0 

Velocity (cm/sec) 



(b) 

Figure 5: Comparison between fuzzy and mathematical model for data 
set II: (a) force-time history, (b) force- velocity. 





Figure 7: Comparison between fuzzy (black) and mathematical 
models (blue) for data set III. 



Chirp displacement- pulse voltage 





Figure 8: Comparison between fuzzy (black) and mathematical 
models (blue) for data set IV. 
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Chirp rib placemen M in iisuUIu 1 1 ullage 




x(k) 



v(k) 



Time (see ) 

Figure 9: Comparison between fuzzy (black) and mathematical 





A quick way to do input selection for inverse model 
is proposed by Roger Jang [35] based on the 
assumption that ANF1S model with the smallest root 
mean square error (RMSE) after one epoch of 
training has a great potential of achieving a lower 
RMSE when give more epochs of training. In our 
study, we tried five models with a different 
combination of four inputs and train them with a 
single pass of the least square method (Table IV). 
Figure 10 shows that the smallest training error is 
achieved when the inputs are force, displacement, 
velocity at the current time step, and voltage from 
previous time step to produce output voltage at the 
current time step. Then the selected inputs are 
trained using the hybrid learning rule to tune the 
membership functions as well (see Figure 11 -a). 

The inverse model is used to predict the 
requirement voltage, which input to the MR damper 
model to produce the desired force as illustrated in 
application phase (Figure 11-b). 



model, a variety of investigation data are applied. 
The investigation data has been collected for five 
seconds and sampling at 2000 Hz. The input is the 
displacement of 2.5 Hz sinusoidal excitation with an 
amplitude 2 cm. 



= ts 

QC 

a 

^ I 
> 

iki 

■ < 

0.01 

£ 0.005 
-s 
> 

iU 0 

GW 

£ 

- 0.01 



■ AN FIS volluge 

■ Target Voltage 



2 3 

Time (sec) 



0.3 



8 0.28 
2 

et 0.26 













i 


f- — I 










► • 1 


I 





Typc(l) Type(2) Type(3) Type(4) Type{5) 

Figure 10: Input selection for inverse fuzzy MR model. 



1 2 3 4 5 

Time (sec) 

Figure 12: Comparison between predicted and target voltage and the 
root mean square errors for constant 1.5V. 

The target voltage is varying with time from 0.5 to 
1.5 V for three tests and constant with 1.5 V for 
fourth test. Figures ( 1 2)-(l 5) compare the predicted 
voltage produced by the inverse fuzzy model with 
the target voltage and present the root mean 
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square error (RMSE) between them. Figure (12) 
and(13) compare the predicted and target when 
the voltage is a constant value of 1.5 V and 
sinusoidal (1+0.5 sin 7it). It is seen that the RMSE is 
about zero and the predicted voltage is match the 
target with high -efficiency. 

Figure (14) and (15) compare the predicted and 
target when the voltage is step and sawtooth with 
an amplitude of 1.5 V. It is shown that the 
predicted voltage is coincided with the target 
except at the peaks with RMSE 0.035 and 0.02 
respectively. 





Time (sec) 



Figure 13: Comparison between predicted and target voltage and the 
root mean square errors for (1+0.5 sin n t) V. 



s 1.5 ■ 



03 r- 



0.04 



C 0.02 



£ 

- 0.02 



0.2 



0.4 ( 1.6 

Time (sec) 



- AN FIS voltage 
■ Target Valin gc 



03 



-<uw 



0.4 0.6 

Time (sec) 



03 



Figure 14: Comparison between predicted and target voltage and the 
root mean square errors for step voltage. 





0 1 2 3 4 5 6 



Time (sec) 



Figure 15: Comparison between predicted and target voltage and the 
root mean square errors for sawtooth with amplitude 2 V. 



5. Conclusion 

The magnetorheological damper is a versatile 
device, which can be used in many applications. The 
identification of MR damper is definitely difficult 
due to its nonlinear dynamics. In this paper, 
ANFIS technique is applied to build the forward 
and inverse fuzzy MR damper 
models. The dynamic behavior of the damper is 
described by the membership function. Data used 
for training the model is generated from numerical 
simulation of nonlinear differential equations. The 
same data are used to train the inverse model, 
which predicts the voltage required to produce the 
desired force. Validation results show that the 
proposed forward fuzzy model can predict the 
damping force accurately. Furthermore, the 
inverse dynamic model can generate the target 
voltage when the MR damper is used in any control 
system. 
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Abstract 

Multilevel power conversion technology is a very 
rapidly growing area of power electronics with good 
potential for further development. Today, it is hard to 
connect a single power semiconductor switch directly 
to medium voltage grids. For these reasons, a new 
family of multilevel inverters has emerged as the 
solution for working with higher voltage levels. 

A Carrier Based Pulse Width Modulation (CBPWM) 
Scheme for multilevel inverters is proposed. The 
proposed PWM scheme generates the inverter leg 
switching times from the sampled reference phase 
voltage amplitudes and centers the switching times 
for the middle vectors, in a sampling interval, as in the 
case of conventional space vector PWM (SVPWM). The 
SVPWM scheme, presented for multilevel inverters, 
can also work in the over modulation range, using 
only the sampled amplitudes of reference phase 
voltages. The sinusoidal PWM technique for various 
modulation indices, Carrier Based PWM technique for 
the cascaded three-level and five-level inverter fed 
induction motor are analyzed and discussed. The 
performance of induction motor is studied in various 
aspects for various modulating techniques and the 
results are compared and tabulated. An induction 
motor model is developed using Simulink and it is fed 
from a three phase three phase five-level cascaded 
inverter. The above said sampled amplitudes of 
reference phase voltages based SVPWM technique is 
implemented using MATLAB/SIMULINK. 

Key words Cascaded Multilevel Inverter, PWM 
Techniques, Carrier Based PWM, Induction Motor, and 
Total Harmonic Distortion. 



1. Introduction 

The general structure of the multilevel converter is 
to synthesize a near sinusoidal voltage from several 
levels of dc levels of dc voltages, typically obtained 
from capacitor voltage sources. As the number of 
levels increases, the synthesized output waveform has 
more steps, which produce a staircase wave that 
approaches a desired waveform [4] [5]. Also, as more 
steps are added to the waveform, the harmonic 
distortion of the output wave decreases, approaching 
zero as a number of levels increases. As the number of 
levels increases, the voltage that can be spanned by 
summing multiple voltage levels also increase. The 
output voltage during the positive half-cycles can be 
found from 

m 
n = 1 

Where SF n is the switching or control function of nth 
node it takes a value of 0 or 1. Generally, the capacitor 
terminal voltages Ei, E 2 ... all have the same value E m . 
Thus, the peak output voltage is v a0 (peak) = (m-1) 
Em=V dc . To generate an output voltage with both 
positive and negative values, the circuit topology [3] 
has another switch to produce the negative part v ob so 
that. 

The multilevel inverters can be classified into three 
types [1, 2]. 

1. Diode-clamped multilevel inverter. 

2. Flying-capacitors multilevel inverter. 

3. Cascade multilevel inverter. 

Among the above specified three topologies Cascade 
Multilevel inverter has the features as follows: 
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For real power conversions from ac to dc and then dc 
to ac, the cascaded inverters needs separate dc 
sources. The structure of separate dc sources is well 
suited for various renewable energy sources such as 
fuel cell, photovoltaic, and biomass. Connecting dc 
sources between two converters in a back-to-back 
fashion is not possible because a short circuit can be 
introduced when two back-to-back converters are not 
switching synchronously [6] [7]. 

The advantages of a cascade multilevel H-bridge 
inverter with s separate dc sources per phase is as 
follows [1, 2] 

• The series structure allows a scalable, modularized 
circuit layout and packaging since each bridge has 
the same structure 

• Requires the least number of components 
considering there are no extra clamping diodes or 
voltage balancing capacitors[8] 

• Switching redundancy for inner voltage levels is 
possible because the phase voltage output is the 
sum of each bridge's. 

Table 1 compares the component requirements per 
phase leg among the three multilevel converters. All 
devices are assumed to have the same voltage rating, 
but necessarily the same current rating. 



Converter 

type 


Diode- 

clamp 


Flying- 

capacitors 


Cascaded 

inverters 


Main 

switching 

devices 


(m-1) x 
2 


(m-1) x 2 


(m-1) x 2 


Main diodes 


(m-1) x 
2 


(m-1) x 2 


(m-1) x 2 


Clamping 

diodes 


(m-1) x 
(m-2) 


0 


0 


Dc bus 
capacitors 


(m-1) 


(m-1) 


(m-l)/2 


Balancing 

capacitors 


0 


(m-1) x 
(m-2) / 2 


0 



Table 1 Comparisons of component requirement per leg of three 
multilevel inverters [2] 



The two most widely used PWM schemes for 
multilevel inverters are the Modified Reference with 
triangular carrier PWM technique [11] and Modified 
Carrier with Sine reference PWM technique. In case of 
first the reference waves can be apart from the sine 
wave i.e. like trapezoidal, step, third harmonic injected 
or space vector PWM etc., and in case of next the 
carrier wave can be apart from the regular triangular 
wave i.e triangle with different phase dispositions and 
u-type carriers with phase dispositions. 



In this paper, the second technique i.e Modified 
Carrier with Sine reference PWM technique and 
compared with the new proposed modulation 
technique called as Modified Reference with Modified 
Carrier modulation technique. In the proposed case 
offset voltage injected reference is compared with the 
u-type carrier signals with different phase dispositions. 
This obtained gate pulses are applied to a Cascaded 
Five level inverter when it feeds to an induction motor. 
This proposed modulation technique called as 
Modified SVPWM, because the reference wave in the 
case of SVPWM looks similar with the proposed one, 
but in performance wise, implementation wise the 
proposed one has no sector identification with vector 
based technique and cumbersome calculations. This 
paper has carried out the simulation results for 
Sinusoidal PWM with different carrier waves and 
Modified Reference with Modified Carrier PWM 
techniques when implemented to a Cascaded five 
level inverter and feed a induction motor. The 
comparative results are taken at output voltage of 
inverter and performance tested in terms of THD. 

2. Modulation techniques 

Sinusoidal PWM 

Pulse width modulation control is most widely used 
method of controlling the modulation depth of 
inverters, including the multilevel family. A significant 
amount of research has been published on the various 
ways of Implementing PWM control. The focus here is 
on carrier-based sinusoidal PWM schemes for 
controlling a cascaded multilevel inverter. 

The SPWM schemes are more flexible and simpler to 
implement, but the maximum peak of the 
fundamental component in the output voltage is 
limited to 50% of the DC link voltage [10], and the 
extension of the SPWM schemes into the over- 
modulation range is difficult. The SPWM technique, 
when applied to multilevel inverters, uses a number of 
level-shifted carrier waves to compare with reference 
phase voltage signals. 

SPWM is the most widely used PWM technique for 
inverters of two and more levels. The basic principle of 
the SPWM used in two level inverters is explained 
here. The reference signal (V r ) which generally is 
sinusoidal is compared with the high frequency 
triangle wave (V c ) of constant amplitude. When V r > 
Vc, the PWM output will be high (state +1), and output 
will be low (state -1) for V r < V c . 

In general, the modulation index is defined as the ratio 
of the magnitude of the reference signal to that of the 
magnitude of the carrier signal. The modulation index 
is m = V r / Vc, Where V r is the magnitude of the 
sinusoidal reference signal waveform and V c is the 
magnitude of the triangular carrier signal. By varying 
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V r and keeping V c constant, the modulation index can 
be changed. By changing the modulation index, the 
amplitude of the fundamental component of the 
output can be varied. For three phase inverters, the 
same carrier signal V c is used for all the three phases 
and three reference signals which are phase displaced 
by 2tt/3 radians are used for each of the phases. 

The parameters of the modulation process are : The 
ratio P = (jo c /u) m , The ratio M - A m and the angle 0 of 
displacement existing between the sinusoidal 
reference and the first positive triangular carrier 
signals 

We also have other degrees of freedom. In fact, the 
phase displacement between two contiguous 
triangular carriers is free. We consider only three very 
simple dispositions that seem the most interesting: 

•All the carriers are phase disposition(PD) 

•All the carriers above the zero value reference are in 
phase among them but in opposition with those 
below (POD) 

•All the carriers are alternatively in opposition 
Disposition (APOD) 

In this paper both triangular and u-type carriers are 
analyzed. The cascaded five-level inverter with R-load 
is simulated and the results are explained below: 

Phase Disposition (PD) 

In this method carriers are the same in frequency, 
amplitude and phases, but they are just different in DC 
offset to occupy contiguous bands as shown in Figure 
1. The carriers are in phase across all the bands. For 
this technique, significant harmonic energy is 
concentrated at the carrier frequency, but since it is 
co-phase component, it does not appear in the line- 
to-line voltage. 



Figure. 1 shows the generation of gate pulses with 
triangular carriers of cascaded five-level inverter with 
R-load. In this four triangular carriers are selected on 
the basis of the formula m-1, i.e 5-1=4. Figure. 2 
shows the generation of gate pulses with u-type 
carrier 

Phase Opposition Disposition (POD) 

Carriers in this method are the same in frequency and 
amplitude but they are again different in phase and 
DC offset as the case of APOD. But in this method 
carriers above the reference zero point are out of 
phase with those below that by 180° as shown in 
Figure.3 




Figure. 3: Generation of gate pulses with triangular carriers 





Tiiiie(ms)-^ 



Figure. 1: Generation of gate pulses with triangular carriers 




Alternative Phase Opposition Disposition (APOD) 

In this method carriers have the same frequency and 
the same amplitude but they are different in their DC 
offset and phases as shown in Figure. 5. In this method 
carriers are phase shifted by 180°, so this method uses 
two degrees of freedom of carriers namely their DC- 
offset and phases. 
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Figure. 6: Generation of gate pulses with U-type carriers 



3. Modified SVPWM 



In SVPWM scheme gives a more fundamental voltage 
and better harmonic performance compared to the 
SPWM schemes [12]. The maximum peak of the 
fundamental component in the output voltage 
obtained with space vector modulation is 15% greater 
than with the sine-triangle modulation scheme. But 
the conventional SVPWM scheme requires sector 
identification and look-up tables to determine the 
timings for various switching vectors of the inverter, in 
all the sectors. This makes the implementation of the 
SVPWM scheme quite complicated the over 
modulation range, extending up to six-step operation. 
It has been shown that, for two-level inverters, a 
SVPWM-like performance can be obtained with a 
SPWM scheme by adding a common mode voltage of 
suitable magnitude, to the sinusoidal reference phase 
voltage. A simplified method, to determine the correct 
offset times for centering the time durations of the 
middle inverter vectors, in a sampling interval, is 
presented, for the two-level inverter. The inverter leg 
switching times are calculated directly from the 
sampled amplitudes of the reference three-phase 
voltages with considerable reduction in the 

computation time [13]. 

To obtain the maximum possible peak amplitude of 
the fundamental phase voltage, in linear modulation, a 
common mode voltage, Voffsetl, is added to the 
reference phase voltages , where the magnitude of 
Voffsetl is given by[13] 






- (V + V . ) 

V max min / 



offset l 



( 1 ) 



In (1), V max is the maximum magnitude of the three 
sampled reference phase voltages, while V^n is the 
minimum magnitude of the three sampled reference 
phase voltages, in a sampling interval. The addition of 
the common mode voltage, Voffsetl, results in the 
active inverter switching vectors being centered in a 
sampling interval, making the SPWM technique 
equivalent to the modified reference PWM technique . 
Equation (1) is based on the fact that, in a sampling 
interval, the reference phase which has lowest 
magnitude (termed the min-phase) crosses the 
triangular carrier first, and causes the first transition in 
the inverter switching state. While the reference phase, 



which has the maximum magnitude (termed the max- 
phase), crosses the carrier last and causes the last 
switching transition in the inverter switching states in 
a two-level modified reference PWM scheme. Thus the 
switching periods of the active vectors can be 
determined from the (max-phase and min-phase) 
sampled reference phase voltage amplitudes in a two- 
level inverter scheme. The SPWM technique, for 
multilevel inverters, involves comparing the reference 
phase voltage signals with a number of symmetrical 
level-shifted carrier waves for PWM generation. 
Because of the level-shifted multi carriers, the first 
crossing (termed the first-cross) of the reference 
phase voltage cannot always be the min-phase. 
Similarly, the last crossing (termed the third-cross) of 
the reference phase voltage cannot always be the 
max-phase. The idea behind the proposed scheme is 
to determine the sampled reference phase, from the 
three sampled reference phases, which crosses the 
triangular first (first-cross) and the reference phase 
which crosses the triangular carrier last (third-cross). 
Once the first-cross phase and third-cross phase are 
identified, the principles of offset calculation of 
equation (1), for the two-level inverter, can easily be 
adapted for the multilevel modified reference PWM 
generation scheme. The proposed modified reference 
PWM technique presents a simple way to determine 
the time instants at which the three reference phases 
cross the triangular carriers. These time instants are 
sorted to find the offset voltage to be added to the 
reference phase voltages for modified reference PWM 
generation for multilevel inverters for the entire linear 
modulation range[14], so that the middle inverter 
switching vectors are centered (during a sampling 
interval), as in the case of the conventional two-level 
modified reference PWM scheme. 

For this analysis we have developed two types of 
carrier based techniques: 

• Modified reference with triangular carriers 

• Modified reference with modified carriers 

As like discussed in SPWM technique the carrier based 
techniques are [11, 13]: 

1. Phase Disposition (PD) 

2. Phase Opposition Disposition (POD) 

3. Alternate Phase Opposition Disposition (APOD) 

Phase Disposition 

In this method carriers are the same in frequency, 
amplitude and phases, but they are just different in DC 
offset to occupy contiguous bands as shown in Figure 
7. The carriers are in phase across all the bands. For 
this technique, significant harmonic energy is 
concentrated at the carrier frequency, but since it is 
co-phase component, it does not appear in the line- 
to-line voltage. 



14 



ACSE Journal, ISSN: 1687-481 1 , Volume 10, Issue 2, ICGST LLC, Delaware, USA, December 2010 




Time(ms)-^ 

Figure. 7: Generation of gate pulses with triangular carriers 







Figure. 10: Output line-line voltage and harmonic spectrum with u- 
type carriers 




Figure. 7 shows the generation of gate pulses with 
triangular carriers of cascaded five-level inverter with 
R-load. In this the reference is modified and four 
triangular carriers are selected on the basis of the 
formula m-1, i.e. 5-1=4. Figure. 3.16 shows the 
generation of gate pulses with u-type carriers. 

Figure. 8 is the line-line output voltage of cascaded 
five-level inverter with R-load and harmonic spectrum. 
For this the triangular carriers are used. The harmonic 
spectrum of cascaded five-level inverter with R-load is 
shown in Figure. 9. It shows that the total harmonic 
distortion is 3.68%. 

Figure. 10 is the line-line output voltage of cascaded 
five-level inverter with R-load and harmonic spectrum. 
For this the u-type carriers are used. The harmonic 
spectrum of cascaded five-level inverter with R-load is 
shown in Figure. 10. It shows that the total harmonic 
distortion is 4.68%. 



400 
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-400 



FFT window: 4 of 25 cycles ilf selected signal 
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Figure. 9: Output line-line voltage and harmonic spectrum with 
triangular carriers 



Phase Opposition Disposition (POD) 

Carriers in this method are the same in frequency and 
amplitude but they are again different in phase and 
DC offset as the case of APOD. But in this method 
carriers above the reference zero point are out of 
phase with those below that by 180° as shown in 
Figure.ll. 




Figure. 11: Generation of gate pulses with triangular carriers 




Figure. 12: Generation of gate pulses with u-type carriers 



Alternative Phase Opposition Disposition (APOD) 

In this method carriers have the same frequency and 
the same amplitude but they are different in their DC 
offset and phases as shown in Figure.13. In this 
method carriers are phase shifted by 180°, so this 
method uses two degrees of freedom of carriers 
namely their DC-offset and phases. 




Tiine(ms)-^ 

Figure. 13: Generation of gate pulses with triangular carriers 
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Among the discussed techniques, PD technique has 
less harmonic distortion on line voltages. As it shown 
PD technique puts the harmonic energy directly into a 
common mode carrier component so that the 
harmonics are cancelled in line voltages. 

4. Simulation results comparisons 

In this section the cascaded multilevel inverter 
simulation results can be shown with R-load by using 
the modulation techniques of SPWM and Modified 
SVPWM. The comparison of simulation results are 
shown in Table 2 & 3. 





Sinusoidal Reference(SPWM) 


Modula 


triangular carrier 


u-type carrier 


PWM 


tion 


PWM 












index 


PD 


POD 


APO 


UP 


UPO 


UAP 








D 


D 


D 


OD 


M = 1 


5.08 


5.87 


12.3 


5.9 


8.81 


12.7 








8 


7 




1 


II 

o 

bo 


5.26 


6.80 


12.9 


6.3 


9.89 


12.8 


66 






0 


0 




8 


II 

o 

bo 


7.36 


9.29 


17.4 


7.1 


11.1 


14.6 








0 


3 


0 


5 



Table 2 The %THD comparison of voltage for SPWM for 
Cascaded Five-level Inverter 



Modula 

tion 

index 


Modified Reference(Modified SVPWM) 


Triangular Carrier 
PWM 


u-type Carrier 
PWM 


PD 


POD 


APO 

D 


UP 

D 


UPOD 


UAPO 

D 


M = 1 


3.68 


5.11 


8.57 


4.6 

8 


4.68 


11.41 


M=0.8 

66 


4.45 


5.23 


10.3 

8 


5.5 

6 


6.32 


12.56 


II 

o 

bo 


4.69 


6.39 


12.2 

5 


6.3 

6 


7.46 


13.01 



Table 3 The %THD comparison of voltage for Modified SVPWM for 
Cascaded Five-level Inverter 



5. CCSLIfed induction motor 

In this Section the simulation results of cascaded five- 
level inverter (CC5LI) fed induction motor, line-line 
voltage of harmonic spectrum, and speed torque 
characteristics are analyzed and discussed. The 
modeling of induction motor as follows. The per phase 
equivalent circuit of the induction motor can be used 



for only steady state analysis[15]. The induction motor 
contains time varying mutual inductances in the 
transient state which makes the analysis using per 
phase equivalent circuit and impossible one[16]. The 
analysis using d-q model can take care of the time 
varying quantities in the transient state. 

The modeling equations in flux linkage form are as 
follows: 




dF dr 

dt 



= co b 





F dr ) 



The complete block diagram of cascaded five-level 
inverter fed induction motor drive is shown in 
Figure. 15. ; in this the cascaded five-level inverter[6] is 
connected at input supply of the system. Two 
comparators are used in the firing circuit, one is for 3- 
phase sine-wave reference signal and the other is for 
level shifted carriers. By using of this combination total 
24 pulses is generated, according this 24 switches are 
required to generate these pulses. 



6 -INDEPENDENT 
DC SOURCES 




Figure. 15: Block diagram of cascaded five level inverter fed 
induction motor drive 
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Figure. 16: Line-line output voltage and FFT spectrum for cascaded 



5-level inverter fed 



IM with SPWM 
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Figure.16 Shows the Line-Line output voltage of 
cascaded five-level inverter fed induction motor. The 
harmonic spectrum is also shown with THD 6.09%. The 
PD modulation scheme of SPWM is used for this 
simulation. 
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Figure. 17: Line-Line output voltage and FPL spectrum for cascaded 
5 -level inverter fed IM with modified SVPWM 



Figure.17 Shows the Line-Line output voltage of 
cascaded five-level inverter fed induction motor with 
no-load. The harmonic spectrum is also shown with 
THD 4.14%. The PD modulation scheme of Modified 
SVPWM is used for this simulation 




Time (s)-^ 

Figure 18 : Stator current of cascaded 5-level inverter fed IM with 
step load T L =7 N-m 

The Output current of induction motor with step-load 
T l =7 N-m is shown in Figure.18. It shows that the 
starting current of induction motor is nearly 12A, and 
after motor reaches the rated speed the current is 
nearly l-2Amp, when load is applied the current 
increases to 5-6 Amp. 




inverter fed with step Load=7 N-m 

Figure.19 shows the Speed and torque waveform of an 
induction motor with step-load of 7 N-m. It shows 



that the speed will reaches the rated speed within 
short time and it retains the rated speed of 1440rpm 
continuously. When a load of 7-Nm is applied, speed 
is reduced and retained 1050 rpm. The torque 
characteristics are also shown in Figure. 19. At the 
time of starting it shows that the torque oscillates, 
when the motor reaches its rated speed it comes to 
constant and when the load is applied on the motor 
the load torque again increases. 

6. Conclusion 

The cascaded multilevel inverter has the superior 
advantages over diode-clamped, capacitor-clamped 
multilevel inverter topologies. So in this paper, the 
cascaded multilevel inverter is proposed for 
implementing of multilevel inverter fed induction 
motor. 

A summary of THD for Cascaded multilevel inverter by 
using SPWM and Modified SVPWM modulation 
techniques is presented. Cascaded 5LI are simulated 
for SPWM, modified reference with triangular carriers 
and modified reference with modified carrier (U-type). 
The modeling of induction motor is developed. For 
this various equations are derived and based on that 
equations the Simulink block diagram of induction 
motor is developed. 

The simulation results of cascaded five- level inverter 
fed induction motor with modified SVPWM are 
analyzed with respect to line-line voltage, the THD 
spectrum of line-line voltage and torque speed 
characteristics of induction motor. The %THD of Line- 
Line output voltage of cascaded five-level inverter fed 
IM with load7N-m by using the modified SVPWM 
Technique is 4.14%. 
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Abstract 

A new method for model reduction is proposed. 
The reduction technique is based a recently 
developed methods, namely frequency domain 
balanced truncation within a frequency bound. This 
reduction technique is developed for model 
reduction of very large dynamical systems which 
are frequently appear in analysis and simulation in 
Micro-technology as a result of finite element 
method. The proposed method is of interest for 
practical model order reduction because in this 
context it shows to keep the accuracy of the 
approximation as high as possible without 
sacrificing the computational efficiency. The 
method is implemented for simulation and analysis 
of micro-electromehcanical optical filter. Numerical 
results show the accuracy and efficiency 
enhancement of the method in comparison to 
previous techniques in this context. 

Keywords: hybrid method, Dynamical systems, 
Frequency-Domain Grammian, Model Reduction, 
Micro-technology. 

1. Introduction 

Over the past two decades model reduction has 
become an ubiquitous tool in a variety of 
application areas and, accordingly, a research focus 
for many mathematicians and engineers. The ever- 
increasing need for accurate mathematical 
modeling of physical as well as artificial processes 
for simulation and control leads to models of high 
complexity. This problem demands efficient 
computational prototyping tools to replace such 



complex models by an approximate simpler models, 
which are capable of capturing dynamical behavior 
and preserving essential properties of the complex 
one, either the complexity appears as high order 
describing dynamical system or complex nonlinear 
structure. Due to this fact model reduction 
methods have become increasingly popular over 
the last two decades. Such methods are designed 
to extract a reduced order state space model that 
adequately describes the behavior of the system in 
question. 

In the field of simulation and analysis in micro and 
nanotechnology frequently we have to deal with 
large dynamical systems which are mostly the 
results of the tools for FEM(finite element method). 
On the other hand these results are in descriptor 
structure and usually singular. This problem is 
motivated the researchers in this context to study 
model reduction of singular systems. In this paper 
an accuracy and efficiency enhanced method for 
model order reduction is presented for model 
reduction of very large singular system. This 
method is a two-step framework. In the first step 
we apply one of the numerically efficient reduction 
techniques to overcome very large orders problem. 
Second step in the framework is frequency-domain 
balanced reduction within frequency bound. This 
framework provides more efficient and accurate 
results comparing to the previous counterparts. 

The paper is organized as follows. In section 2, we 
present some system theoretic and mathematical 
background briefly. Section 3 consists of presenting 
our method in detail. In section 4 the proposed 
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methods applied to the model in microtechnology. 
The conclusion is presented in section 5. 

2. Mathematical and System Theoretic 
Background 

In this section some notions and notations in 
context of model reduction in general and model 
reduction of singular systems in particular are 
presented in brief. More details can be found in [1]- 
[4]. 

Def. 1. Dynamical singular systems: 



This kind of transformation enable us to get non- 
singular system out of a singular structure ( by 
expressing x 2 in terms of X x and eliminating x 2 ). 

3. Frequency-Domain Hybrid Reduction 
Method 

In this section we first present frequency domain 
balanced truncation method and the related 
techniques and then the frequency-domain hybrid 
reduction method for reduction of very large 
dynamical systems. 



Ex = Ax +Bu, y =Cx (1) 

and 

Ex = Ax +Bu, y - Cx (2) 

are called restricted equivalent systems if there 

exist non-singular P and Q such that: 

x=Px, QEP = E, QAP =A, QB=B , CP =C 



Def. 2. Dynamical singular system (1) is regular if 

there exist a eC such that: 

det (aE + A ) ^ 0 or det (aE - A ) ^ 0 

Following we recall a theorem from [4] which is 
very important in reduction of dynamical singular 
systems. 



Theorem 1. Every regular singular system has a 
restricted equivalent dynamical system which is 
described by the following equations: 



x i — A^x i A^x 2 B 


(3) 


0 = A 2\X i + A 22 X 2 + B 


(4) 


y =C 1 x 1 +C 2 x 2 


(5) 


In order to get this structure we can 


apply SVD on 



E: 

U t EV =diag(L, 0) 

and then by coordinate transformation: 



We have: 



x =V 



X , 






A, 

_ Ai 




[C, C 2 ] = cv 1 




TU t B 

o" 

/ 



3.1 Model order reduction within 
frequency bound 

The Consider the following n' h order state-space 
model representation of an asymptotic stable LTI 
system: 



fx') 




r A B ' 








D y 


Uy 



( 6 ) 



The problem consists in approximating the system 



with r th order state-space model: 



(A r ,B r ,C r ,D r ) 



(7) 



where r < n . 

There are huge amount of algorithms are available 
in the literature for model reduction, which some of 
them are summarized in Figure.l. 



One commonly used and globally more accurate 
approach is the so-called Balanced Model 
Reduction first introduced by Moore [7][19][20]. In 
this method, the system is transformed to a basis 
where the states which are difficult to reach are 
simultaneously difficult to observe. Then, the 
reduced model is obtained simply by truncating the 
states which have this property. Two other closely 
related model reduction techniques are Hankel 
Norm Approximation [11] and the Singular 
Perturbation method [12]. When applied to stable 
systems, all of these approaches are guarantees to 
preserve stability and provide bounds on the 
approximation error. Because of being operational 
the system within frequency band and outside that 
it is not important to have an accurate 
approximation the accuracy can be improved with 
applying balanced model reduction in the specified 
frequency band [8],[9][14-18],[22], 

Controllability and observability Gramians in terms 
of w over a frequency band [w 1? w 2 \ are defined as 
[8],[9][14-18]: 
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Figure!. Algorithms for model/controller reduction. 

W, := — P ( Ijw - A )" 1 BB* i-Ijw - A* f' dw 

2ft Jw 2 

W of : = P { ~ ] i w ~ ) _1 C f C(Ijw - A)~ l dw 

] 2ft Jw 2 



Those are the solutions for the following Lyapunov 
equations: 



AW cf +W cf A* = -(BB*F* + FBB *) 



A*W of +W of A = -(C*CF + F*C*C) 



Where F is defined by: 

F = P 2 (Ijw- A)' 1 dw 

J vv, 

With an appropriate similarity transformation T and 
change of the basis, system realization in (7) can be 
transformed to a new balanced realization, so that 
the Gramians are equal and diagonal (in decreasing 
diagonal elements): 

W cf = W of =diag( (T l ,(T 1 ,...,(T n ) 



Here we have two important theorems that give us 
a physical interpretation for the reduction 
procedure [8],[9][14-18]: 



Theorem 2:The frequency-domain controllability 
Gramian represents the energy flow of the system 
through each state variable within the frequency 
range [w v w 2 ] . 

It means that if the unit white Gaussian noise test 
input signal u(t) , and state vector x(t)of the system 
defined as follows: 

, ,2 fl w, -< I w| -< w 9 

u(jw)[u(jw)Y = \u(jw)\ =<^ 

[0 elsewhere 

x(t) = e A, Bu(t ) jwl - A) 1 Bu(jw) = X(jw ) 

The energy of the system through controllability 
Gramian is as follows: 

E x = x(t)x (t)dt = 

— f 2 (jwl - Ay 1 Bu(jw)u* (jw)B* (-jwl - A*)~ l dw 

2ft Jw i 

= W cf [w v w 2 ] 



Theorem 3: The frequency-domain observability 
Gramian represents the energy flow of the system 
through each state variable within the frequency 
range [w, , w 2 ] . 

Consider a unit injected test signal x 0 , where 

* fl w^<w<w 2 
x 0 x 0 = < 

[0 elsewhere 

Define output 

y(t) = Ce A x 0 <=> C (jwl - A ) -1 * 0 = Y (>v) 



The energy E of the system through observability 
Gramian is obtained by: 

E y ={ o y(t)y(t)dt = 

— I* 2 (-jwl - A*)~ l C*C(jwI - A)~ l dw 
In 



Now, x 0 being a white noise test signal, the result 
follows:^ =W cf [w v w 2 ] 

From the above theorems, it is understood that for 
having a good approximation we should only 
truncate the states that are related to the lowest 
singular values. 



3.2 Mixed method for model order 
reduction 



In order to have more flexibility and more accuracy, 
the singular perturbation method and the 
proposed frequency balance method are combined. 
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The mixed method is based on the fact that every 
linear time-invariant system can be decomposed 
into the sum of two subsystems (slow and fast) in 
frequency domain. 

System (7) can be partitioned into the slow and fast 
modes as: 







(a a b ^ 

A 1 A 2 






x 2 (t) 


= 


Al A2 A 




x 2 (t) 


K y(t)) 




K C\ C 2 D y 


f 


v u(t) , 



To consider model reduction at a particular 
frequency levelG, we apply singular perturbation 
and replace x 2 (t) by crx 2 (t) so that G can be 
chosen to be any frequency within [w p w 2 ] . 

Setting s = ja and applying the method the 
reduced order model is obtained as: 



(a) 




Al + A 2 ~ A2) Al 


B 




A A2 ( ~ A2 ) A 


C 




Q (j&In-r ~ A2 ) Al 


, D . 




V ^ ) + Q(7 <T A-r “ A2) A > 



'f 



The Figure (2) shows the flowchart for the above 
method: 




Figure 2. The mixed method Flowchart. 



The idea of I/O weighting balance method is first 
proposed by Enns. Figure(3) shows the schematic 
of the I/O weight applied to a system with transfer 
function G(s ): 




Figure 3. The schematic of I/O weights. 



To combine the I/O weighting balance method and 
the proposed method in Figure (1) for an input or 
output weight W- or W Q and the system with 
realization: 



‘A 


B ” 




"A B.~ 




‘A b~ 


0 


O 


,W:= 


l l 


,G(s ) := 




C 0 


D 

0 J 


C, D t _ 


C D 



The overall transfer function of the structure in 
Figure (3) should be obtained: 



G(s)^(s) = 



'A 


B~ 


[4 


Bi] 


C 


D_ 


1% 


D >\ 



A Bq 


BDi 


0 Ai 


% 


C DQ 


DI\ 



~A 


S' 


C 


D_ 



In the similar procedure we could obtain the overall 
transfer function for the output weights. 

The overall transfer function obtained above is 
used as an input for the method proposed in 
Figure (2). The accuracy and performance is 
increased when we use the proposed mixed 
method and the proposed technique preserves and 
improves the advantages of the singular 
perturbation, frequency balanced and the I/O 
weighting reduction method. 



3.3 Practical consideration 

One should know how to deal with repeated or 
almost equal singular values in frequency-balanced 
realization when applying the frequency balanced 
truncation method. Like the usual balanced 
truncation. If we have two equal singular values and 
we only truncate the state which is related to one 
of them the stability of the original system may not 
be preserved. A large amount of the computational 
load of every Gramian based reduction method is 
imposed by finding the Gramians. In the proposed 
model reduction method we can easily obtain the 
Gramians by using the definitions and 

approximating the integral to summation. 
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3.4 Benchmark examples for frequency 
domain balanced reduction 

In this section of the paper we apply the presented 
method to some benchmark examples to illustrate 
it more. The examples are chosen from the large 
scale "real word" systems reflecting current 
problems in applications [13]. 

Building Model 

One of the important benchmark examples in the 
context of the linear time invariant systems is the 
building model of the Los Angeles University 
Hospital with 8 floors each having 3 degrees of 
freedom, namely displacement in x and y directions, 
and rotation. Hence we have 24 variables with the 
second order differential equation describing the 
dynamics of the system. The system is single input 
single output and the order of the system in the 
state-space form is 48. We apply usual balanced 
method and the proposed method to this model 
and reduce the order of it to 4. The infinity norm 
error of the method is shown in Figure (4). The error 
is lower than usual balanced method with more 
efficiency in computations. 



Approximation Errors 




Figure 4. The approximation error in usual balanced method(solid 
line) and the error in the proposed method in the frequency band 
[0.1,4] (dash dotted). 

CD Player Model 

We apply the proposed method to another 
benchmark. This system describes the dynamics 
between the lens actuator and the radial arm 
position of a portable CD player. The model has 
120 states with 2 inputs and 2 outputs. The systems 
is reduced to 8 th order by the balanced truncation 
method and reduced to 4 th order by the proposed 
method. The approximation error is shown in figure 
(5).AIthough the order of the reduced system by 
the proposed method is lower than the order of the 



reduced system by balanced truncation technique, 
the reduced model by proposed balance method is 
more accurate than the reduced system by usual 
balanced truncation in some frequency ranges. 



Approximation Errors 




Frequency (rad/sec) 

Figure 5. The approximation error in usual balanced 
method(solid line) and the error of the proposed method in the 
frequency band (dash dotted). 

Atmospheric Storm Track 

The model of an atmospheric storm track is a large 
scale single input single output model with 548 
states. The system is reduced to 4 th order model and 
the results are shown in figure (6). The infinity norm 
error of the proposed method is lower than the 
infinity error of the balanced truncation reduction 
method. 



Approximation Errors 




Figure 6. The approximation error in usual balanced method(solid 
line) and the error in the proposed method in the frequency band 
[0.9,20] (dash dotted). 

3.5 Frequency-Domain Hybrid 
Reduction Method 

At this point we are ready to develop new model 
reduction technique for reduction of large- 
dynamical system. First we use theorem 1 to the 
general structure: 
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Ex(t)=Ax(t) + Bu(t), y(t)=Cx(t ), t> 0 
and we achieve non-singular system : 

x | = (A^ — j + (B ^ — A 12 A 2 2 ^ 2 )^ 

y = (C 1 — C 2 ^ 22 ^ 21 )^ 1 "I - ( — 1 ^ 2 ^ 22 ^ 2 ^ 

which is ready for reduction. In the first step we 
apply one of the methods from moment matching 
based reduction[6]. The second step would be 
applying frequency-domain balanced truncation 
within frequency domain which is described in 
detail in [10]. 

4. Benchmark from micro-technology: 
Tuneable Optical Filter 

In this section we apply our method to a 
benchmark example in microtechnology and we 
compare the proposed method with the hybrid 
method which is based on the method in [3]. 

Model order reduction approach is successfully 
used to considerably reduce the computational 
time and resources. Mathematical development of 
MOR is an active area of research, which grows 
from the 

reduction of linear ordinary differential equation 
systems (ODEs) towards the reduction of 
parametrized and nonlinear differential algebraic 
equations (DAEs) and partial differential algebraic 
equations (PDAEs). The implementation aspects of 
model order reduction are advancing as well. 
Practical MOR has developed from academic 
prototyping environments to several strong tools 
that can be easily used as an extension of the 
commercial simulators like e. g. ANSYS[10]. 

There are several electro-thermal MEMS models 
that can show the application of model reduction 
in Micro-Technology. We have chosen tunable 
optical filter to apply our method [21]. 

The DFG project AFON (funded under grant ZA 
276/2-1) aimed at the development of an optical 
filter, which is tunable by thermal means. The thin- 
film filter is configured as a membrane in order to 
improve thermal isolation. Fabrication is based on 
silicon technology. Wavelength tuning is achieved 
through thermal modulation of resonator optical 
thickness, using metal resistor deposited onto the 
membrane. The devices features low power 
consumption, high tuning speed and excellent 
optical performance. 

Figure(8) shows the structure of one of the 
industrial benchmark examples in microtechnology. 
Optical tunable filter is a micro-electromechanicall 
system which is tunable by thermal means. The 
benchmark contains a simplified thermal model of 
a filter device. It helps designers to consider 



important thermal issues, such as what electrical 
power should be applied in order to reach the 
critical temperature at the membrane or 
homogeneous temperature distribution over the 
membrane. 

The original model is the heat transfer partial 
differential equation.The output that we are 
interested in is the temperature in the center which 
is shown with 1 in the Figure (10). More details 
regarding the technical issues can be found in[10]. 
Figure.7 shows the procedure for using model 
reduction in simulation and analysis of 
Microsystems. 




Figure 7. General Procedure for analysis and simulation for 
MEMS 

Figure. 9 shows the procedure for using model 
reduction in simulation and analysis of 
Microsystems in our case study. 




Figure 8. Tunable optical filter 



If we apply FEM(finite element method), we end up 
with a SISO dynamical system of order 1668 which 
really hard to simulate and analysis. According to 
our proposed framework we first reduce this 
dynamical system to a dynamical system of order 
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200 and in the next step we apply frequency- 
domain balanced truncation. 

Singular Values 





Figure 9. Procedure for analysis and simulation for MEMS in our 
case study 




Figure 10. Cross Sectional view with outputs. 



Figure 11. Singular values: original model(solid), reduced 
method by the proposed method(dotted) and method based on 
[3](dash-dotted) 



Singular values 




Figure 12. Singular values: original model(solid), reduced 
method by the proposed method(dotted) and method based on 
[3](dash-dotted) 



Step Response 




Figure 13. Step responses: original model(solid), reduced 
method by the proposed method(dotted) and method based on 
[3](dash-dotted) 



Singular value of the reduced model and original 
systems are depicted in Figure(ll) and Figure(12). 
In Figure (11) the frequency domain hybrid method 
has reduced the dynamical system to a system of 
order 7 th and the method based on [3] to 20 and in 
figure (12) frequency domain hybrid method has 
reduced the dynamical system to a system of order 
7 th and the method based on [3] to 15. It is clear 
from the Figures that the proposed framework for 
reduction provides better results and can 
approximate the dynamical system into lower 
orders while preserving accuracy of the 

approximation. Figure(13) shows the step 

responses of the original and reduced models 
which confirms the accuracy of the approximation. 
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5. Conclusion 

In this paper a framework for model order 
reduction of large dynamical systems with 
descriptor or singular structure is proposed. This 
method is a two-step framework. In the first step 
we apply one of the numerically efficient reduction 
techniques to overcome very large orders problem. 
Second step in the framework is frequency-domain 
balanced reduction within frequency bound. This 
framework provides more efficient and accurate 
results comparing to the previous counterparts. 
This reduction technique is developed for model 
reduction of large dynamical singular systems 
which are frequently appear in analysis and 
simulation in Micro-technology as a result of finite 
element method. Numerical results show the 
accuracy and efficiency enhancement of the 
method in comparison to previous techniques in 
this context. 
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Abstract 

This paper presents an analytical approach for the 
analysis of the current ripple in space vector based 
PWM controlled induction motor drives. This paper 
considers seven different types of switching 
sequences. By using the instantaneous ripple 
voltages, the current ripple expressions are derived 
for all the sequences. To reduce the complexity 
involved in the conventional space vector approach, 
the proposed sequences have been carried out 
using the concept of imaginary switching times. 
Finally, to reduce the current ripple, hybrid PWM 
(HPWM) algorithm has been developed which uses 
a sequence that results in the lowest rms current 
ripple over the given sampling time interval. Thus, 
the proposed HPWM algorithm reduces the current 
ripple at all modulation indices. To validate the 
proposed HPWM algorithm, numerical simulation 
studies have been carried out and compared with 
the conventional space vector PWM (CSVPWM) 
algorithm. 

Keywords: current ripple , hybrid PWM, induction 
motor, SVPWM. 

1. Nomenclature 

Van ,V b n,V cn : Instantaneous phase voltages 

Tan > Tbn > T cn : Imaginary switching times 

7j , 7 2 : Active vectors switching times 

T z : Zero voltage vectors switching time 

^max’^min ^mid '• Maximum, minimum and middle 

values of imaginary switching times 

V rip : Ripple voltage vector 



V re f : Reference voltage vector 

a : Angle of the reference voltage vector 
i rip : Current ripple vector 

i rms : Rms value of current 
My. Modulation index 
cr : Leakage coefficient of induction motor 
L s : Stator inductance in H 

2. Introduction 

Recently, the voltage source inverters are becoming 
popular to feed variable voltage, variable frequency 
voltages to the three-phase induction motor drives. 
A detailed survey on the various PWM algorithms 
for voltage source inverters is given in [1]. Among 
the various PWM algorithms, SVPWM algorithm 
offers many advantages and hence it is becoming 
popular [2]. Moreover, among the various 
approaches, the triangular comparison (TC) 
approach and space vector (SV) approach are 
popular. However, the SV approach offers so many 
degrees of freedom for the implementation [3]. 
However, the conventional SVPWM (CSVPWM) 
algorithm uses equal division of zero state time 
duration. By utilizing the freedom of active state 
division, nowadays many PWM algorithms are 
derived [3]. By using only one zero state, many bus- 
clamping PWM algorithms can be generated [4]. But 
the existing bus-clamping PWM algorithms, which 
are given in [4] did not involve in the active state 
division. Similar to the zero state time division, the 
active state time division also can be performed in 
the SV approach. Recently, these are attracting the 
researchers and a detailed study about these 
algorithms are given [3], [5]-[7], [9], [11]. The double 
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switching sequences (involving active state division) 
give superior performance at higher modulation 
indices. To reduce the steady state ripples at all 
modulation indices hybrid PWM (HPWM) algorithms 
are developed by using the notion of stator flux 
ripple [5]-[7], [9], [11]. 

However, the conventional SV approach requires the 
calculation of angle and sector information in every 
sampling time interval and hence increases the 
complexity. To reduce the complexity involved in 
the conventional SV approach, a novel approach has 
developed by utilizing the concept of offset time in 

[10] . To reduce the burden involved in conventional 
HPWM algorithms, novel HPWM algorithm using 
the concept of imaginary switching times is given in 

[11] . However, the proposed HPWM algorithms in 
[5]-[7], [11] uses the notion of stator flux ripple to 
develop the HPWM algorithm. To reduce the 
current ripple, the analysis has been carried out 
using the notion of the current ripple in [8]-[9]. 

This paper presents a novel HPWM algorithm for 
reduced current ripple which uses the notion of 
stator current ripple and uses the concept of 
imaginary switching times [10]-[11]. 

The remainder of the paper is organized as follows: 
Section (2) focuses on the generation of various 
switching sequences using the concept of imaginary 
switching times. Section (3) emphasizes on the 
analysis of stator current ripple. Section (4) focuses 
on the development of HPWM algorithms and 
section (5) gives the simulation results and 
discussions and conclusions are given in section (6). 

3. Switching Sequences 

The conventional SV approach uses reference frame 
transformation and angle calculation. Hence, it 
increases the complexity involved in the algorithm. 
But, the proposed approach uses the instantaneous 
reference phase voltages for the generation of 
switching sequences. For the explanation purpose, 
the concept of imaginary switching times will be 
introduced. The imaginary switching times are 
proportional to the instantaneous phase voltages 
and can be defined as follows: 

T an = TT- V an ’ T bn = 77 ^ V bn ’ T cn = T7~ V cn (!) 

v dc V dc V dc 

Yaw Vb n an( ^ V cn are the instantaneous reference 
phase voltages and T an , T bn and T cn are the 
corresponding imaginary switching times. As these 



times are proportional to the instantaneous voltages, 
these times could be negative where the voltages 
are negative. Hence, these times are defined as 
imaginary switching times. The active vector 
switching times Ti and T 2 , if the reference voltage 
vector falls in sector-1 may be expressed as follows 
[ 11 ]: 

T l = T an ~ T bn * T 2 = T bn ~ T cn ( 2 ) 

In the conventional SVPWM algorithm, when the 
reference voltage vector falls in the first sector, the 
imaginary switching time which is proportional to 
the a-phase ( T an ) has a maximum value, the 
imaginary switching time which is proportional to 
the c-phase ( T cn ) has a minimum value and the 
imaginary switching time which is proportional to 
the b-phase (T bn ) is neither minimum nor maximum 
value. Thus, in general to calculate the active vector 
switching times, the maximum (r max ), middle ( T mid ) 
and minimum (r min ) values of imaginary switching 



times are calculated in every sampling time. Then 
the active vector switching times Ti and T 2 may be 
expressed as 



7 \ 7" max T m f^ , T 2 T mid T mm (3) 

The maximum and minimum imaginary switching 
times and active vector switching times for all six 
sectors are given in Table 1. 

The zero voltage vectors switching time is calculated 
as 



T z ~T S -T\- T 2 (4) 

To generate the various switching sequences, the 
zero state time will be shared between two zero 
states as xT z for V 0 and (1 -x)T z for V 7 respectively. 
Thus, by using (3) and (4), the active vector times 
and zero vector times can be calculated. If x=0.5, 
then the conventional SVPWM algorithm is obtained. 
When x=0, any one of the phases is clamped to 
positive dc bus for 120 degrees over a fundamental 
cycle. When x=l, any one of the phases is clamped 
to negative dc bus for 120 degrees over a 
fundamental cycle. Then by utilizing the space 
vector concept and zero and active vector division 
concept, various space vector based PWM 
algorithms can be generated. The existing bus- 
clamping (BCPWM) algorithm use only one zero 
state and the advanced BCPWM (ABCPWM) 
sequences use only one zero state and active state 
division. The sequences of all the possible PWM 
algorithms in six sectors are listed in Table-2. 
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Table 1: Various times in all six sectors 





Sector-I 


Sector-II 


Sector-Ill 


Sector- 

IV 


Sector-V 


Sector- 

VI 


T 

1 max 


T 

± an 


T bn 


T bn 


T 

± cn 


T 

± cn 


T 

± an 


^min 


T 

* cn 


Ten 


Tan 


T 

± an 


T 

1 bn 


T 

1 bn 


7 ] 


T T 

1 an 1 bn 


T T 

1 bn 1 an 


T T 

1 bn L cn 


T T 

1 cn 1 bn 


T -T 

*cn ± an 


T -T 

± an *cn 


TV 


T T 

1 bn 1 cn 


T -T 

* an A cn 


T -T 

*cn ± an 


T T 

1 bn 1 an 


T T 

1 an 1 bn 


T T 

1 cn 1 bn 



Table 2: Sequences of various possible PWM algorithms 



Sector 


CSVPWM 


BCPWM algorithms 


ABCPWM algorithms 


I 


0127-7210 


012-210 


721-127 


0121-1210 


7212-2127 


1012-2101 


2721- 

1272 


II 


0327-7230 


032-230 


723-327 


0323-3230 


7232-2327 


3032-2303 


2723- 

3272 


III 


0347-7430 


034-430 


743-347 


0343-3430 


7434-4347 


3034-4303 


4743- 

3474 


IV 


0547-7450 


054-450 


745-547 


0545-5450 


7454-4547 


5054-4505 


4745- 

5474 


V 


0567-7650 


056-650 


765-567 


0565-5650 


7656-6567 


5056-6505 


6765- 

5676 


VI 


0167-7610 


016-610 


761-167 


0161-1610 


7616-6167 


1016-6101 


6761- 

1676 



4. Analysis of the Current Ripple 

The space vector approach generates the reference 
voltage vector in an average manner. But, when the 
required voltage vector is applied to the motor, the 
current is flowing in the stator winding of the 
induction motor. The actual value of the stator 
voltage differs from the applied voltage. Hence, 
there is always an instantaneous error voltage vector. 
The error voltage vector is defined as given in (5) 

y rip =v k -v ref ( 5 ) 

where 'k' is the k th voltage vector and it takes any 
value from '0' to 7'. Then, it is possible to define the 
current ripple vector as defined in (6) 




The voltage ripple vectors and the variation of the 
current ripple by the application of corresponding 
voltage vectors are shown in figure 1 for all the 
possible PWM algorithms. Also, from this figure, it 
can be observed that the current ripple has zero 
mean value in a sampling time period. Then the 
current ripple vectors are given by as follows: 



hi ~ 



Vl ~ Vref_ 

crL s 



T , 



(7) 



ir 2 - 



Vl-Vref 

cfL s 

aL s z 



t 2 



( 8 ) 

(9) 



The current ripple vectors can be represented in 
terms of imaginary switching times as given in (10)- 
(12). 

'Vl = id + jiq\ ( 10 ) 

'V 2 = id + jiq2 (11) 



l rz = ~J 



T- 2M:V, 



i v dc 



oL c 



- = ~Ji 



qz 



( 12 ) 



where MAs the modulation index and defined as 



M: =- 



7TV \ 



ref 



Wdc 
„ ( 



l d 



l q\ 



oL v 



7 lV f 



dc 



3y[3 M t T 



ciL, 



l q2 



s J 

2Vy c ^(7i+0-5r 2 ) IVj'Mj 
9M t T s n 

r 2V dc 7T(0.5 T l+ T 2 ) 2 V dc Mf 



CTL r, 



(13) 



(14) 



(15) 



(16) 



9M t T s 7i 

Then the current ripple vectors can be resolved 
along the d- and q- axes, which are the 
synchronously revolving reference frame axes as 
shown in figure 2, from which it can be observed 
that the application of a zero state results in a 
variation of the q-axis component of the current 
ripple and the application of any active state results 
in variation of the both the d-axis and q-axis current 
components. 
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V 2 





(b) 



V 2 




(c) 





Figure 1: Error voltage vectors along with the current ripple 
variations (a) 0127, 012, 721 sequences (b) 0121 (c) 7212 (d) 
1012 (e) 2721 




(a) 




(b) 




- 0-5 id 



(c) 
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.2 



rms$)\2\ 3 



i ilz. 
l qz T 

1 S 

iqz + i* 



iqz^qz + Q-^ql) + i^qz + Y J 



2T ' 



P iqz + 5iql Y kqz + ^ql F* ^ql + (p* ^ql Y J 

(o-^lP^r+lO.Sy) 2 ^^) 



(d) 




Figure 2: the q-and d-axes components of the current ripple 
vectors for the possible sequences, (a) 0127, 012, 721 sequences 
(b) 0121 (c) 7212 (d) 1012 (e) 2721 

The q-axis, d-axis and the total rms current ripple 
over a sampling time interval (T s ) can be evaluated 



as follows: 




1 T s 

.2 1 p.2 t 

l q,rms ~ T J l q,rip ^ 
0 


(17) 


1 T s 

• 2 1 p .2 7 

i A = — \i A . dt 

a, rms 7 ^ J a, rip 

0 


(18) 


•2 -2 -2 

l rms =l q,rms +l d,rms 


(19) 



Finally, the rms current ripple expressions can be as 
given in (20) - (24). 



•2 xT z 

XI r , — ^ + 

«° T s 



2 . 2 
xiqz + ( iq Z + iq\) + iqz^qz + iqD 



5-+ 

JT, 



(iqz + *q0 0 - x )iqz(iqz + iq0 + (Q x ^qz) 

(d - X )iqz) 2 ^-+id- Tl+T2) 



T s 



T s 



T 2 

— + 



( 20 ) 



•2 



1 



rm$12\T 3 



l qz T 



,• 2 . 

l qz 



+ hjkqz + d^/2 ) + {iqz + ®--*q2 F 



( 21 ) 



2 T s 



{iqz + ^q2^ {iqz + ®-^q2p-^q2 + (d-5^r2p 1 1 



T 2 ,(nc-.\2^+^) 



K 2 p^-+(o.%) : 

zi s 



Ts 



( 22 ) 



,2 J 

W012T 3 



+ p0-5/^l) 2 +05i q i[o.5i q i+i q ^+ip.5i q i+i q ^ J- 
+p0.5/^j +/<^p +(0-5/^i +^)+(^l +^zF 

+fel +'J 2 y +i d — - +(°-5y ) 2 f 1 

' V A if 



(23) 



,- 2 -I 

rm&llT 3 



+[(o.5/^2p + 03 ^ 2 ( 0 . 5^2 + z ^)+(o5^2 +^z) 
+p0.5/^2 +^zP +(o3^2 +iqzfaq2 +i q) + ^q2 +i qzj^ 



(24) 

The number of switchings of the CSVPWM and 
ABCPWM sequences in a sampling time interval is 3 
and whereas for the BCPWM algorithms is 2. Hence, 
to get the constant average switching frequency of 
the inverter, a sampling time interval is taken as 
T s = T for the CSVPWM and ABCPWM algorithms, 
while T s =(2r/3) for the BCPWM algorithms. Thus 



the proposed method of analysis is capable of 
estimating d-axis and q-axis stator current ripples 
individually. Total harmonic distortion (THD) of the 
current waveform is equally affected by the d-axis 
and q-axis stator current ripples. On the other hand, 
torque pulsation mainly depends on q-axis current 
ripple, and is practically independent of d-axis stator 
current ripple. Thus, the reduction in q-axis current 
ripple signifies the reduction of THD and torque 
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pulsation, while the reduction in d-axis current 
ripple implies reduction in THD only. By using (20)- 
(23), the rms current ripple trajectories over a sixty 
degrees period for every PWM sequence at different 
modulation indices are plotted as shown in figure 3. 




(a) 




(b) 




(c) 




(d) 

Figure 3: current ripple trajectories over a sixty degrees period 
(for one sector) (a) Mf0.4, (b) Mf0.6, (c) Mf 0.8 and (d) 
MfO. 905 (A=0127, B=012, C=721, D=0121, E=7212, F=1012 and 
G= 2721) 

5. Development of HPWM Algorithms 

From figure 3, it can be observed that replacing of 
a by (60 ° -a) in the current ripple expression of 
CSVPWM sequence does not change their values. 
Thus for a given sampling time and M { , 0127 
sequence produces equal mean square current 
ripple at spatial angles a and (60 0 -a). Whereas 
for the other sequences the mean square current 
ripple at a is equal to the mean square stator 
current ripple of their respective complimentary 
sequences at (60 ° -a) . For a given values of 
M i and T s , these can be given as follows: 



irms, 012 (^) — * rms, 121 (^0 — Oc) 


(25) 


irms, 0121 (^) — ^rms,1212 (^0 ~ &) 


(26) 


i rms, 1012 (^) — ^ rms, 2721 ~ 


(27) 



By comparing the sequences 0127, 012, 721, 0121, 
7212, 1012 and 2721 with respective to each other, 
the zones of superior performance of each 
sequence can be identified. The development of 
HPWM techniques for reduced current ripple 
involves determination of superior performance for 
every sequence. The zone of superior performance 
for a given sequence is the spatial zone within a 
sector where the given sequence results in less 
mean square current ripple than other sequences 
considered. Thus, the proposed HPWM algorithm is 
a 7-zone PWM algorithm, which uses six bus- 
clamping PWM algorithms along with the CSVPWM 
algorithm for reduced steady state current ripple. 

The flowchart of proposed HPWM algorithm is as 
shown in figure 4. 
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6. Simulation Results and Discussion 

To verify the proposed HPWM algorithm, numerical 
simulation studies have been carried out on a v/f 
controlled induction motor drive by using 
Matlab/Simulink. For the simulation, the average 
switching frequency of the inverter is taken as 5 kHz. 
The induction motor used in this case study is a 4 
KW, 400V, 1470 rpm, 4-pole, 50 Hz, 3-phase 
induction motor having the following parameters: 
R s = 1.57Q, R r = 1.210, L s = 0.17H, L r = 0.17H, 
L m = 0.165 H and J = 0.089 Kg.m 2 . The steady state 
currents along with their harmonic spectra are given 
in Figure. 5 - Figure. 8 for SVPWM and proposed 
HPWM algorithms based v/f controlled induction 
motor drives at different supply frequencies. From 
the simulation results, it can be observed that the 
proposed HPWM algorithm reduces the total 
harmonic distortion (THD) in the line current when 
compared with the CSVPWM algorithm. Moreover, 
from the line current spectra corresponding to 
CSVPWM algorithm, it can be observed that the 
harmonic components are around frequencies that 
are integral multiples of the switching frequency (5 
kHz). On the other hand, the proposed HPWM 
algorithm gives spread spectra. Hence, the 
proposed HPWM algorithm also reduces the 
acoustical noise of the induction motor with a 



uniform sampling frequency. Thus, the proposed 
HPWM algorithm reduces the steady state current 
ripple and acoustical noise of the motor when 
compared with the CSVPWM algorithm. 




0 5000 10000 15000 20000 

Frequency (Hz) 

(b) 



Figure 5: Steady state results with SVPWM algorithm at a supply 
frequency of 35 Hz, (a) Steady state current (b) current spectra 




(a) 



ACSE Journal, ISSN: 1687-481 1 , Volume 10, Issue 2, ICGST LLC, Delaware, USA, December 2010 




Frequency (Hz) 
(b) 



Figure 6: Steady state results with SVPWM algorithm at a supply 
frequency of 45 Hz, (a) Steady state current (b) current spectra 
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Figure 7: Steady state results with HPWM algorithm at a supply 
frequency of 35 Hz, (a) Steady state current (b) current spectra 
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0 5000 10000 15000 20000 



Frequency (Hz) 

(b) 

Figure 8: Steady state results with HPWM algorithm at a supply 
frequency of 45 Hz, (a) Steady state current (b) current spectra 

7. Conclusions 

A simplified HPWM algorithm is presented in this 
paper for reduced current ripple. This paper has 
analyzed the instantaneous current ripple and 
derived the expressions for the variation of current 
ripple for the possible PWM algorithms. Finally, 
HPWM algorithm has developed in this paper. From 
the simulation results, it can be observed that the 
proposed HPWM algorithm gives less THD when 
compared with the conventional SVPWM algorithm 
over a wide range of operating frequencies 
(modulation indices). Moreover, as the HPWM gives 
spread spectra, it reduces the acoustical noise of the 
induction motor with a uniform sampling frequency. 
Hence, to reduce the acoustical noise, special 
methods with random sampling frequencies are not 
necessary. Thus, the proposed HPWM algorithm 
gives better performance at all modulation indices. 
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Abstract 

In chemical industries parallel cascade controllers 
are used to improve the dynamic performance of a 
control system. Cascade controllers has two 
controllers, primary and secondary controllers in 
cascade. In this paper, a fuzzy based control is 
designed for the primary controller to improve the 
performance of the parallel cascade control. The 
secondary controller is implemented with internal 
model control (IMC) method. The dead time in the 
process is compensated by incorporating Smith 
predictor in the primary loop. The use of fuzzy logic 
controller is more effective than the conventional 
controller. It contributes to achieve quick response 
and high accuracy to variations in input parameter 
fluctuations during operations in systems, and also 
provides robust control performance. The 

performance of the proposed method has been 
evaluated by conducting experiments on two 
different systems. The simulation results revealed 
that the proposed method outperforms other 
recently reported methods in performance 

measures for both transient and steady state 
conditions. 

Keywords: Parallel cascade control, PID control, 

Fuzzy logic control, Dead time, Internal model control 

1. Nomenclature 

G p i, G p2 Primary and secondary processes 

transfer function 



Gdi, G d 2 Transfer functions of disturbances (d) 

Gd, G c2 Primary and secondary controller 

transfer function 

d Disturbance input for primary and 

secondary loop 

ri, r 2 Primary and secondary loop set point 

Yi, Y 2 Primary and secondary process output 

G p2m Secondary process model 

Gf 2 Set point filter 

Gfi Filter for predicted disturbance 

G pm Process model of the outer loop 

without dead time 

0 m Dead time of the model of the outer 

loop 

u Crisp valued controller output 

fi Fuzzy valued controller output 

Pi Input membership function 



2. Introduction 

In conventional single feedback control, the 
corrective action for disturbance does not begin 
until the controlled variable deviates from the set 
point. If the response is unsatisfactory because of 
large uncontrolled load changes, more complex 
control schemes can be considered [1], Cascade 
control, first introduced by Franks and Worley [2], is 
one of the frequently used schemes to deploy the 
availability of additional controllers for the plant to 
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effectively adopt to wide fluctuations in loads. The 
concept of cascade control is to nest one feedback 
loop inside another loop. Cascade control involves 
the use of a primary and secondary controller. In 
process industries, cascade control is widely used to 
reduce the effects of possible disturbances and to 
improve the dynamic performance of the cascade 
loop system. In traditional series cascade control, 
both the manipulated variable and the disturbance 
affect the primary output through the intermediate 
(secondary) output. The cascade control of parallel 
type arises when the manipulated and disturbance 
variables simultaneously affect primary and 
secondary outputs [3]. The overhead composition 
control of a distillation column, cascaded onto the 
control of a tray temperature, is a typical example of 
a parallel cascade control. The reflux flow rate 
(manipulated variable) and the feed flow or 
composition (disturbance) affect both, the purity of 
the overhead product (primary output) and the tray 
temperature (secondary output). The control 
objective is to maintain the overhead composition 
at set point. The output of the composition 
controller resets the set point for the temperature 
controller. By controlling the tray temperature in the 
cascade manner, the variation in the feed can be 
compensated before it disturbs the product 
composition. Cascade control is used to improve the 
dynamic response of the column [4]. In general, 
parallel cascade control is appropriate when the 
secondary loop has a faster dynamic response and 
the rejection of the disturbance in the secondary 
output reduces the steady state output error in the 
primary loop. The parallel cascade control is also 
beneficial when measurements of the primary 
output are sampled infrequently and/or with long 
time delays. 

Though parallel cascade control is widely used in 
industries, it has attracted very little research despite 
its clear benefits. Yu [5] proposed an interaction 
measure to determine whether the parallel cascade 
control is advantageous or detrimental to load 
response. Shen and Yu [6] applied these results to 
the selection of secondary measurement in parallel 
cascade control. Bramilla and Semino [7] introduced 
a nonlinear filter in order to practically separate the 
dynamics of the primary and secondary loop in 



parallel cascade control. They also developed a 
combined structure to avoid interactions between 
primary and secondary loops [8]. Mcavoy and Ye [9] 
discussed nonlinear inferential parallel cascade 
control structure based on both parallel cascade 
control and inferential sensing. Pottman et al. [10] 
developed a parallel control strategy for a biological 
control system. Chen et al. [11] focused on the 
performance assessment of parallel cascade control 
systems using the methodology of univariate 
control loop performance, minimum variance and 
Diophantine decomposition principles. Lee et al. [12] 
proposed an analytical method of Proportional 
Integral Derivative (PID) controller design for parallel 
cascade control system which takes into account the 
interaction between primary and secondary loops. 
Recently, Seshagiri et al. [13] incorporated time 
delay compensator in the existing parallel cascade 
control. Most industrial processes and the primary 
loop in cascade control systems in particular, are 
non-linear, multivariable, distributed complex 
dynamic systems. There are also large delays in 
dynamic channels and cross relationship between 
control parameters. The process dead time and gain 
are non-linear and are not predictable because they 
vary according to process conditions [14]. For 
instance, changes in flow, pressure, temperature, or 
level may cause large variations in process dead 
time and gain. This makes some process control 
loops extremely non-linear. Due to 

interconnections, dead time and non-linearity, the 
control of these systems become difficult [15]. Also, 
the accurate process models are not known. When 
there is a lack of accurate process model, the 
existing controllers designed so far will not be much 
effective. Fuzzy control strategy is fit for systems 
that lacks accurate model. In this paper, fuzzy 
control is designed for the primary controller. It is a 
well known method of implementing nonlinear 
control. Hence the parallel cascade control is 
improved by incorporating fuzzy controller in 
controller design. 

This paper is organized as follows: Next section 
gives brief overview of parallel cascade control. 
Section 4 provides the secondary and primary 
controller design procedure and simulation result 
for the different case studies to illustrate the use of 
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proposed cascade control. Conclusions are provided 
in section 5. 

3. Parallel Cascade Control 

The parallel cascade control structure is shown in 
Figure 1. G p i and G p2 are the transfer functions of the 



primary and secondary processes respectively. Gdi 
and G d2 are the transfer functions of disturbances 
(d) for the primary and secondary loops. G c iand G c2 
are the primary and secondary controllers 
respectively. 




Fig lira 1 , Parallel cascade control system 



In parallel cascade control, the manipulated and 
disturbance variables simultaneously affect primary 
and secondary outputs. The control objective is to 
maintain the primary output Yi at the set point in 
the midst of disturbances. Y 2 is controlled by the 
dynamics of the internal loop. The disturbance 
rejection in the primary loop depends on the 
disturbance rejection in the secondary loop as well 
as its set point tracking. Output of the primary 
controller, G c i resets the set point of the secondary 
loop. The secondary measurement, Y 2 is designated 
to infer the load disturbance. Secondary controller 
G c2 is designed for disturbance rejection purpose, 
while the primary controller is designed to maintain 
good servo response. Hence Parallel cascade control 
structure offers an effective solution to overcome 
disturbance problems. 

The proposed parallel cascade control is shown in 
Figure 2. In parallel cascade control, the dynamics of 
the secondary loop is faster than the primary loop. 
They have negligible or no dead time compared to 
primary process, which has larger dead time. So 
secondary controller, G c2 is designed based on IMC 
method. G p2m is the secondary process model and 
G f2 is the set point filter. The primary process are 
usually non linear with large dead time and they are 
difficult to model. Non-linear boundary conditions, 
time varying characteristics and a complex physical 
environment create many difficulties in designing a 



mathematical model. So, the existing PID or IMC 
based controllers will not be effective. They hardly 
satisfy a high performance criterion with respect to 
both parameter insensitivity and disturbance 
rejection. Also, effectiveness of IMC controller 
depends on the accuracy of the process model. 
Though it gives good servo response, its load 
regulation is poor. On the other hand, fuzzy control 
is one of the useful control techniques for uncertain 
and ill-defined nonlinear systems because available 
qualitative operator and design knowledge can be 
easily implemented. Hence the existing control can 
be improved by incorporating fuzzy controller for 
the primary controller, G c i. Further, the dead time in 
the outer loop is compensated by incorporating 
Smith predictor in the primary loop, which 
compensates the dead time [16]. G pm is the process 
model of the outer loop without dead time. 0 m is the 
dead time of the model. By including the filter G f i 
with time constant T f , for the predicted disturbance, 
the robustness of the delay compensator can be 
increased [17], where 



Gn - 



TfS + 1 



( 1 ) 



The closed loop transfer functions of secondary and 
primary loop for the system in Figure 2 are, 

for secondary loop (assuming G p2 = Gp 2m ) 

Y 2 = G p2m G c2 G f2 r 2 + C 1 “ G p2m G c2 ) G d2 d ( 2 ) 



41 



ACSE Journal, ISSN: 1687-481 1 , Volume 10, Issue 2, ICGST LLC, Delaware, USA, December 2010 



Further, (since G p2 = Gp 2m ), the relation between 
primary output (YJ and secondary set point (r 2 ) is 
obtained as 

Y 1 = G pl G c2 G f2 r 2 



However, for the primary loop, we also have 

y _ Gf 2 G pl G cl G C 2 ^ 

1 l+G c i(G pm (l- Gfxe -0 !^ s ^_|_G pl G c2 Gf! Gf 2 ) 1 

( 1+G pm G ci( 1- G fi e_0m S ))( G di -G pi G c2 G d2) 

— - Q (4) 

1+G c i (G pm (l- G^e -0 ™ s s)+G pl G c2 G^ Gf 2 ) 




From Eq. (3), by assuming perfect model conditions, 
the primary process model is obtained as 

Gpm = G p iG c2 G f2 (5) 

4. Controller Design and Simulation 

Simulation results are obtained for different case 
studies. Experiment is carried out for both slow and 
fast dynamic systems. The performance of the 
designed controller is compared with the recently 
reported methods. 

Case study No.l 

The following process was modeled by Semino and 
Bramilla [8], The process is characterized by the 
same dynamics in G pi and G p2 and in G d i and G d2 
with exception of a dead time which corresponds to 
a delay in measurement of the primary output. 



Gp 2 



- Gd2 - 



3.14e~ 9s 
30s + 1 



(7) 



Here the primary process is dead time dominant 
and the secondary process is purely lag dominant. 
Lee et al. [9] have designed the secondary controller 
for slow dynamics and the primary controller for the 
fast dynamics and found improved results. For the 
proposed method the Internal model controller is 
designed for the secondary controller and fuzzy 
logic controller is designed for the primary 
controller. 



IMC based Controller Design (G c2 ) 

As stated in section I, the secondary controller is 
designed based on IMC principles [12, 18]. The 
inner loop transfer function in the form of a First 
order plus time delay (FOPTD), can be expressed as 



G 



P2 - 



k 2 e ® 2 s 

Cr 2 s + l) 



( 8 ) 



The process and disturbance transfer function 
models are given as 

l.24e -33s 



Where, k 2 is the process gain and t 2 is the process 
time constant. Then the secondary controller using 
IMC principle, is given by 



G p i - G dl 



30s + 1 



( 6 ) 
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_ (t 2 s + 1) 
c2 k 2 (A s + 1) 



(9) 



where, A is the only tuning parameter to be found. 
Smaller the value of A, the better the performance of 
the cascade control system. If a faster response is 
required, A can be chosen as half time delay of the 
inner loop, i.e., A = 0 2 /2. The lead term in the closed 
loop response of the secondary output for set point 
changes may cause an overshoot and large settling 
time which will affect the closed loop system 
performance. To eliminate this undesirable 
overshoot, a set point filter G f2 is used [12]. 



The secondary controller tuning parameter is 
selected as A = 5, so as to achieve the same 
secondary disturbance rejection as that of Seshagiri 
Rao et al.[13]. The secondary controller is obtained 
as 



G c2 - 



(30s + l)(14.6s + l) 



3.1(5s+l) 2 

and the loop filter 

Gf2 — 

rz 14.6 s + 1 

The primary process model is obtained as 



Gpm e 



_q s 0.4e 
D m b — 



( 10 ) 



(ii) 



( 12 ) 



(5 s + l) 2 

Using Skogestad's half rule [19], the above model is 
reduced to FOPTD Model and is obtained as 



G 



pm 



e 



0 46-35-5 S 

(7.5S+1) 



(13) 



Fuzzy Logic Controller Design (G cl ) 

Control algorithms with fuzzy logic controllers (FLC) 
prove to offer better robustness and efficiency in 
the case of more complicated non-linear and time 
varying systems and systems with large dead time in 
comparison with conventional PID controllers. While 
a PID controller allows no flexibility of structure, a 
fuzzy logic controller can be whatever it needs to be 
[20]. The difficulty in obtaining the optimum settling 
time of PID controller is mitigated by using fuzzy 
logic controller which gives the opportunity to 
describe the control action in quantitative form and 
symbolic form. Hence, fuzzy logic controller is 
designed for the primary controller (G c i). It must be 
recalled that the general design of the fuzzy logic 
controller proceeds along the following steps: 



(i) Fuzzification 

(ii) Design of Membership function 

(iii) Defuzzification 



The design process needs a systematic method for 
obtaining the rule base and domain ranges. The 
error (e) and change in error (Ae) are chosen as an 
input to FLC. 
where, 

Error = Setpoint - Measured Variable (14) 



Change in Error = Present Error- Previous Error (15) 
Since the design of FLC is based on "If - Then" Rule- 
based considerations, it must be noted that in the 
present context, the two signals given in (12) and 
(13) constitute the rule-antecedent in the formation 
of rule base [i.e., the "If"-part]; the control output (u) 
represents the contents of the rule-consequent in 
the formation of rule base [i.e., the "Then"-part]. 
Fuzzification 

Fuzzification is the quantization of crisp quantity to 
fuzzy linguistic variables. In many of the quantities 
that are considered crisp and deterministic are 
actually not deterministic at all. Fuzzification is 
related to the vagueness and imprecision in a 
natural language. It is a subjective valuation, which 
transforms a measurement into a valuation of a 
subjective value. Hence, it could be defined as a 
mapping from an observed input space to fuzzy sets 
in certain input universes of discourse. Fuzzification 
plays an important role in dealing with uncertain 
information, which might be objective or subjective 
in nature [21]. 

Design of Membership function 

The ranges of input and output are assigned with 
linguistic variables. These variables transform the 
numerical values of the input of the fuzzy controller 
to fuzzy quantities. These linguistic variables specify 
the quality of the control. Among all membership 
functions, triangular membership function gives 
good results and easy to implement. In this work, 
the triangular membership function is used for 
designing the FLC. The designed fuzzy triangles for 
Error (e), Change in Error (Ae) and the controller 
output (u) are shown in Figure 3. 

The Linguistic variables for error (e) is classified into: 
Negative big (e_ big ); Negative Small (e_ sm ); Zero (e- 
zero ); Positive Small (e +sm ); Positive Big (e +big ) 

The Linguistic variables for Change in Error (Ae) is 
classified into: 

Negative big (Ae. big ); Negative Small (Ae. sm ); Zero 
(Ae zer o); Positive Small (Ae +sm ); Positive Big (Ae +big ) 

The Linguistic variables for FLC output (u) is 
classified into: 

Negative big (u _ big ); Negative Small (u _ sm ); Zero 
(Uzero); Positive Small (u +sm ); Positive Big (u +big ) 

Knowledge Base 

The knowledge base of an FLC is comprised of two 
components, a database and fuzzy control base. The 
concepts associated with a database are used to 
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characterize fuzzy control rules and a fuzzy data 
manipulation in FLC. It should be noted that the 
correct choice of the membership function of a term 
set plays an essential role in the success of the 
application. A lookup table based on discrete 
universes defines the output of a controller for all 
possible combinations of the input signals. A fuzzy 
system is characterized by a set of linguistic 
statements. It is in the form of "IF-THEN" rules. 
Fuzzy conditional statements form the rules or the 
rule set of the FLC. The decision table is given in 
Table 1. 



Table 1: Rule Matrix of the proposed FLC 



Ae(k) 


6-big 


6-sm 


6zero 


6+sm 


6+big 


Ae_big 


U -big 


U -big 


U -big 


U -sm 


Uzero 


A 6 -sm 


U -big 


U -big 


U -sm 


Uzero 


U +sm 


A^zero 


U -big 


U -sm 


Uzero 


U +sm 


U +big 


A£+sm 


U -sm 


Uzero 


U +sm 


U +big 


U +big 


A^+big 


Uzero 


U +sm 


U +big 


U +big 


U +big 




Figure 3: Fuzzy triangles for Error (e), Change in Error (Ae) and 
Controller output (u) for Case study No. 1. 



Defuzzification 

The final stage of fuzzy logic implementation is 
defuzzification. In this process, the fuzzy values are 
converted into acceptable crisp values. Centroid 
defuzzification is used in the majority of the cases to 
develop the expected value for each of these output 
variables [22], 

By the centroid defuzzification method, the crisp 
controller output (u) is obtained by, 



N 

I f. p. 

i ^ 1 1 



u = - 



N 



X 

i = l 



y L ; 



(16) 



where, 

u = Crisp valued controller output 
fi = Fuzzy valued controller output 
Pi = Input membership function 
N = Total number of rules 

From the above formulae, controller output is 
obtained which is given to the secondary process. 
By this way, fuzzy logic controller is designed and 
implemented for the primary controller. 

With the above controller design, the proposed 
method is simulated by giving unit step change in 
the set point at time, t=0s and a negative step input 
of magnitude 2 in the load disturbance at t=350s 
respectively. The closed loop response is shown in 
Figure 4. The closed loop response of the proposed 
method reaches the set point with less rise time and 
settles to the final steady state value with minimum 
overshoot when compared to that of [13]. The 
values obtained are listed in Table 2. Also, it is 
observed that the disturbance rejection of the 
proposed controller is better in terms of fast 
recovery and minimum overshoot. 
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To analyze the robustness of the controller, a 
perturbation of +20% in the primary process time 
delay and +20% of disturbance time delay are given 
to the model. The corresponding closed loop 
response is shown in Figure 5. For obtaining the 
result analysis under transient conditions rise time 
(t r ), peak overshoot (%OS) and settling time (t s ) are 
compared. For obtaining the result analysis under 
steady state conditions integral square error (ISE) 
and integral of absolute error (IAE) are compared. 




Figure 4: Closed loop response for Case study No.l for perfect 
model. Solid - Proposed method, dash - Seshagiri Rao et al. [13]. 




the peak overshoot is reduced by 62% when 
compared to [13]. The ISE and IAE of the proposed 
method are also reduced considerably when 
compared with existing methods. Similar 
improvements were obtained for disturbance 
rejection also. Seshagiri Rao et al. [13] have already 
shown that their method performs better than Lee 
et al. [12]. 



Case study No.2 

Here, an example from Luyben [3] is considered. The 
dead time in the primary process and disturbance is 
considered as 4 times the time constant. The 
resulting process and disturbance transfer function 
models are given as 



G d1 = G dl = 

P 1 ai 20s + 1 



G D 2 — G d2 — 

P z az 10s + 1 



(17) 

(18) 



Here the secondary process dynamics is faster than 
the desired closed loop response. The tuning 
parameter for the secondary controller is selected as 
A = 1, so as to achieve the same secondary 
disturbance rejection as that of [13]. 



The secondary controller is obtained as 

r _ (10s+l)(1.9s+l) 
c2 (s+1) 2 

And the loop filter 

Gf2 — 

IZ 1.9 s + 1 



(19) 



( 20 ) 



The primary process model is obtained as 



G e -9 m s _ (10s+De- 80s 

P m (20s+l)(5 s + l) 2 



( 21 ) 



Figure 5: Closed loop response for Case study No.l for 
perturbation of 20% in the primary process and disturbance 
delay. Solid - Proposed method, dash - Seshagiri Rao et al. [13]. 



Table 2: Performance analysis for Case study No.l 





tr (S) 


%OS 


ts (S) 


ISE 


IAE 


Proposed 

Method 


15.8 


1 


60 


3.78 


17.28 


Seshagiri 
Rao et al 


31.9 


2.6 


90 


12.64 


26 



The values of these different performance indices 
are provided in Table 2. There is 51% improvement 
in rise time, 33 % improvement in settling time and 



Using Skogestad's half rule, the above model is 
reduced to FOPTD Model and is obtained as 



Gpm e 



-0m S _ _JL 



(20.5S+1) 



( 22 ) 



The primary fuzzy controller is designed as 
explained above. The designed fuzzy triangles for 
Error (e), Change in Error (Ae) and the controller 
output (u) are shown in Figure 6. With the above 
controller design, the proposed method is simulated 
by giving unit step change in the set point at time, 
t=0s and a negative step input of magnitude 10 in 
the load disturbance at t=250s respectively. The 
closed loop response is shown in Figure 7. 
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Here too, the closed loop response of the proposed 
method reaches the set point with less rise time and 
settles to the final steady state value with minimum 
overshoot when compared to that if [13]. Also, 
disturbance rejection of the proposed controller is 
better and the system recovers faster and reaches 
the set point with minimum overshoot. Hence it is 
observed that the proposed controller gives better 
response and the control action of the proposed 
method was found to be superior to that of [13]. 
The values of different performance indices are 
provided in Table 3. There was 44% improvement in 
rise time and 3% improvement in settling time. 
There is no peak overshoot when compared to that 
of [13]. The ISE and IAE of the proposed method are 
also very minimum and comparable with existing 
methods. Similar improvements were obtained for 
disturbance rejection also. 

Table 3: Performance analysis for Case study No.2 





tr (S) 


%OS 


ts (S) 


ISE 


IAE 


Proposed 

Method 


16.5 


0 


114 


6.89 


19.53 


Seshagiri Rao 
et. al 


29.4 


8.5 


118 


14.69 


28.91 





Figure 6: Fuzzy triangles for Error (e), Change in Error ((Ae) and 
Controller output (u) for Case study No. 2 




Figure 7: Closed loop response for Case study No.2 for perfect 
model. Solid - Proposed method, dash - Seshagiri Rao et al. 
[13]. 

5. Conclusion 

Fuzzy controller is designed and implemented for 
the primary controller of parallel cascade control 
system. Dead time compensator is incorporated in 
the primary loop. The proposed structure brings 
together the best merits of the fuzzy control and the 
parallel cascade control structure. The performance 
of the proposed method is analyzed for both slow 
and fast dynamic systems and it is shown by the 
case studies, that the fuzzy based parallel cascade 
control gives a better performance for both set 
point tracking and load disturbance rejection. The 
result is a stable control system with fast reaction 
time and the proposed method outperforms the 
recently reported methods in performance 

measures for both transient and steady state 
conditions. 
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Abstract 

A model for the steady state, geometrically non- 
linear, periodic forced vibration of fully clamped thin 
rectangular composite plates under harmonic 
external excitation is presented. The equations of 
motion have been derived using spectral analysis, 
Lagrange's equations, and the harmonic balance 
method (HBM). This forced case has been examined 
using a single and multimode approach. A 
multidimensional Duffing equation in the case of 
fully clamped rectangular laminated composite 
plates has been obtained. The behaviour of the 
symmetrically laminated plates on the responses to 
forced vibration under harmonic distributed force. 
Results obtained are compared with some solutions 
based on hierarchical finite element method. It is 
found that accurate results may be produced using 
the present model. 

Keywords: Clamped thin rectangular composite 
plates Eexternal excitation , Ftarmonic balance 
method ' Multimode approach 

1. Introduction 

With the increasing use of laminated composite 
materials in structures, engineers and researchers 
are faced with challenging problems associated with 
the understanding of the behaviour of these new 
materials. The correct and effective use of them 
requires more complex analysis tools in order to 
predict accurately their response to external 
loadings. As most modern laminated plates are 



made of high performance materials, they are 
usually used in severe working conditions resulting 
in large vibration amplitudes response. This means 
that the effect of geometric non-linearity needs to 
be included. A considerable amount of research 
work has been carried out on geometrically non- 
linear behaviour of laminated composite plates, 
particularly free vibration of thin plates [1 to 10]. 
However, attention has been also devoted to the 
forced vibration of composite plates, as may be 
seen by the amount of published works on the 
subject. In reference [11 to 13] the finite element 
method has been extended to investigate large 
amplitude forced vibrations of rectangular 
composite plates. Based on the hierarchical finite 
element method (HFEM), a forced vibration analysis 
of symmetrically laminated plates has been studied 
in reference [14]. The loads considered are harmonic 
acoustic plane waves impinging on the plate surface 
in normal direction and at grazing incidence. Only 
linear vibrations have been analysed. Recently, a 
model for the steady state, geometrically non-linear, 
periodic forced vibration of both isotropic and 
composite laminated rectangular plates under 
harmonic external excitation has been presented in 
[15]. The model has been derived by applying the 
principle of virtual work and the HFEM. The 
equations obtained have been transformed into the 
frequency domain by the HBM and have been 
solved by a continuation method. In the above 
reference, the convergence properties of the model 
have been discussed by applying it to isotropic and 
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to composite laminated plates. The stability study 
and analysis of multi-modal response of the model 
have been investigated in reference [16]. 

In a previous work [17] a model for the steady state, 
geometrically non-linear, periodic forced vibration 
of fully clamped thin rectangular composite plates, 
under harmonic external excitation was presented. 
The model has been derived using spectral analysis, 
Lagrange's equations, and the harmonic balance 
method (HBM). A multidimensional Duffing 
equation in the case of fully clamped rectangular 
laminated composite plates has been obtained. It 
was assumed for simplicity only one mode, this set 
reduces to the Duffing equation, very well known in 
one mode analyses of non-linear systems having 
cubic non-linearties. The frequency response curves 
have been obtained at the plate centre, for various 
levels of loading. It appeared that the method works 
well, since excellent agreement is found between 
the result of the present model, and those published 
in the literature in the centre of the plate. In this 
paper, the model developed in [17] is extended to 
investigate the multimode response. The response 
of the plates to external harmonic excitations is 
analysed and compared with some published results. 
Fully clamped boundaries are considered because 
they are adequate to model many real panel-type 
situations, such as aircraft wing panels [18]. Results 
are given for single and multimode approach. 

In this paper we will rely on the Lagrange equations 
to study the vibrations of a rectangular plate and 
built a composite material of known characteristics. 
Then we will use Harmonic Balance Method coupled 
with the routine NS01A to develop the method 
mentioned above. And finally we will compare 
results with other studies that already exist 

2. Theory 

Consider a thin plate of thickness H having a length 
a and a width b, and composed of thin orthotropic 
layers, layered along length and width, that is x and 
y axis, alternately. The origin of the Cartesian co- 
ordinate system is located in one corner of the mid- 
plane with the z axis being perpendicular to this 
plane as shown in Figure 1. Assuming that the plate 
undergoes large amplitudes, the Von-Karman type 
strain-displacement relationships as given below are 
used: 



{e} = {e 0 } + z{k} + R°} 



( 1 ) 
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( 2 ) 



(3) 



(4) 



U, V and W are the displacements of the plate mid- 
plane, in the x, y and z directions respectively. 

For the laminated plate having n layers shown in 
Figure 2, the stress in the k th layer can be expressed 
in terms of the laminate middle surface strains and 
curvatures as: 



K}=[ei-M (5) 

In which <7y T - xy J and terms of 

matrix bJ can be obtained by the relationships 
given in [19]. 

The force and moment resultants for the k th layer 
are defined as follows: 




( 6 ) 



( a ,yy)zd Z 

h-i ( 7 ) 

Where h k is the distance from the mid-plane to the 
surface of the k th layer. 

The in-plane forces and bending moments in a plate 
are given by: 



in which {c 0 }, {k} and {: X °} are given by: 
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A, B and D are the symmetric matrices given by (9). 
[B] = 0 for symmetrically laminated plates [20]. 

H 

2 

( A ij ’ B ij ’ D ij)= \QiP(}-’Z’Z 2 )* 

-H 

2 (9) 



Q {k) 

Where the lJ are the reduced stiffness 
coefficients of the k th layer in the plate co-ordinates. 
The expression for the bending strain energy V b , 
axial strain energy V a and kinetic energy T are given 
by [21]: 
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( 11 ) 



( 12 ) 



Where S is the plate surface [0,a] x [0,b] and dS is 
the elementary surface dx dy. 

In the above expressions, the assumption of 
neglecting the inplane displacements U and V in the 
energy expressions has been made as for the fully 
clamped rectangular isotropic plates analysis 
considered in [22,23]. The range of validity of this 
assumption has been extensively discussed in the 
light of the experimental and numerical results 
obtained for the frequency amplitude dependence 
and the bending stress estimates obtained at large 
vibrations amplitudes. The results obtained via this 
assumption were compared with the previous ones 
based on various methods such as the finite 
element method, the method based on Berger's 
approximation, the ultraspherical polynomial 
method and the elliptic function method. It was 
found that the percentage error in the non-linear 
frequency estimates based on this assumption, for 
amplitudes up to 1.5 times thickness, does not 



exceed 1.3 %. Also, in the experimental investigation 
of the non-linear behaviour of fully clamped 
rectangular plates at large vibration amplitudes 
presented in reference [24], it was found that the 
rate of increase in bending stresses estimates, 
obtained from measured data, was in very good 
agreement with that obtained from the theory, in 
which the assumption of zero in-plane 
displacements was made. The above references 
concerned with the isotropic case. It has led to the 
conclusion that this assumption allows a great 
simplification in the modelling and a great reduction 
in the computation time when calculating the non- 
linear mode shapes, the associated non-linear 
frequencies, and the non-linear bending stress 
patterns, for a reasonable range of vibrations 
amplitudes and plates aspect ratios. In the present 
work, concerned with symmetrically laminated 
plates, the validity of this assumption has also been 
examined, via comparison of the non-linear 
frequency estimates it leads to, with those obtained 
by the HFEM, for two symmetrically laminated plates, 
with different lay-up and thicknesses. 

If the time and space functions are supposed to be 
separated and harmonic motion is assumed, and the 
generalised parameterisation and usual summation 
convention defined in [25], the transverse 
displacement can be written as: 

W(x,y,t) = q i (t)w i (x,y) 



where Wj (x,y) are the basic functions. Substituting W 
in the expressions (10-12) for V a , V b and T, the 
discretisation of the kinetic energy and the total 
strain energy V leads to: 



n =y < h<lj k ij 

rji 1 

T = ~q i q j co mij 



(14) 

(15) 

(16) 



where the terms mj, ky and by k | are given as in 
reference [6 to 9, 25], and indices i, j, k and I are 
summed over 1 to n 



2.1 Non-linear Free vibration 
The numerical model developed in [22,26,27], and 
used in [6 to 10], by applying Hamilton's principle to 
the free response problem, can also be derived 
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using Lagrange's equations and the Harmonic 
Balance Method (HBM). 

For a conservative system, having n degrees of 
freedom corresponding to n parameters q r (t), r = 1 
to n, Lagrange's equations can be written, if no 
forcing term is considered, as [28]: 



d 8T | 8T 8V _ Q 
dt dq r dq r dq r 



r = 1 to n 



(17) 



Substituting equations (14-16) in (17) leads to the 
following set of non-linear differential equations, 
similar to that found in [6 to 10,22,29,30], and 
considered as a multidimensional form of the very 
well known Duffing equation: 

q^r +q i K+2q i q j q k b ijkr =0 r=1 tQ n (lg) 
which can be written in matrix form as: 
[Mfel+Kfel + 2[5({«})fe} = {0} 

(19' 



where [M], [K], [B] and {q} are the mass, the linear 
rigidity, the non-linear rigidity tensors and the 
column vector of generalised parameters 



{qf =[q x ,q 2 . 




respectively. 



in which Wj(x,y) is the ith mode of the plate 
considered. 

The non-dimensional axial and bending stiffnesses 

4* D* 
and U *: 
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(23, 24) 



K U TTl 

The non-dimensional tensors ij , ijkl and ij can 
be obtained as: 



aH 5 E * 
k » = — [ k “ 

m kj = pH 3 abm*j 

and: 



co 1 = 



hee 

pb A 



-co 




aH 5 E * 
—^ h m 



(25-27) 



(28) 



Inserting equations (22-28) into equation (21), the 
following equations is obtained: 

|r ]- a- 1 [m • |a} + 1 [zr ({a})](a} = {o} 

2 (29) 

which may be written also, using tensor notation as: 



By assuming that harmonic motion takes place: 

q, (0 = a t cos {cat) (2Q) 

Inserting equation (20) into equation (19), and 
applying the HBM one arrives at the following 
system of equations: 

2([/f]-ffl 2 [M])(A}+3[s({A})]{A} = {0} (21) 



in which {A} is the column vector of the basic 
function contribution coefficients 

=[a x ,a 2 , ,a n ] 

In order to formulate the vibration problem of fully 
clamped-clamped rectangular laminated plates 
corresponding to large deflections in non- 
dimensional terms, we put, as in reference [6,25]: 



w ( . (x, y) = Hw* — , yl = Hw* (x\/) 



a b 



( 22 ) 



~(o i2 a,i>i +a i k* r +-a i a J a k b* kr =0 

1 r = 1 to n (30) 

For the non-linear problem (29), as the modulus of 
(A] is also necessary to determine the physical 
amplitude of vibration, we have (n + 1) unknowns 
which are the n coefficients of {A} and co. A further 
equation has to be added to equations (29) in order 
to complete the formulation. As no dissipation is 
considered here, such an equation can be obtained 
by applying the principle of conservation of energy, 
which can be written, as in reference [25]: 

Vmax = T max (31) 

This leads to the following equation: 



CO = 



a^jk*! + a i a j a k a l b i 



ijkl 



apjTTiy 



(32) 

To obtain the first non-linear mode shape of the 
plate considered, the contribution of the first basic 
function is first fixed and the other basic function 
contributions are calculated via numerical solution 
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of the remaining (n-1) non-linear algebraic 
equations. 



2.2 Non-linear forced vibration 

The model developed above is extended in this 
subsection to the steady-state periodic forced 
response case. The excitations considered will be of 
the form F(x,y,t) = F coscot. Numerical, experimental 
and analytical investigations confirm that plates 
vibrating with displacement amplitudes of the order 
of their thickness due to harmonic excitation forces, 
usually experience periodic motion [29]. In this work, 
only the first harmonic motions is examined, via use 
of the HBM. 



Consider a plate subjected to a transverse 
concentrated force F c applied at the point (x 0 , y 0 ), or 
to a transverse distributed harmonic uniform force 
F d , distributed over the range S, where S is the 
surface of the plate or a part of it, F c and F d are given 
by: 



F c (x, y, t) = F c cos(cot)S(x - x 0 )S(y - y 0 ) (33) 

F d (x,y,t) = F d cos(ctf) (34) 

8 is the Dirac function. The corresponding 

can be 



generalised forces ^ and 



F-(t) 



written as: 



F- (t) = F c cos (cot)w i (x 0 ,y 0 ) = f. cos(cot) 

F d (t) = F d cos(ryf)J w i (x,y)dxdy 



(36) 



= f d cos(ryf) 



f- 



The non-dimensional generalised forces Ji and 
ff 

J ' corresponding to a concentrated force at (x 0 , 
y 0 ), and a uniformly distributed force are given by: 



f- =F C 



aEH 



-w*(Wo) 



l F 



u pp 

fi d ' =F d -j^\yvt(x*,y*)dx*dy* 



(37) 



(38) 



Adding now the forcing term {F(t)} to the right hand 
side of equation (19) leads to: 



[M]{q}+ K\q)+2[B{{q\%q) = {Fit)} 



(39) 



Inserting equation (20) into (39) and applying the 
HBM one obtain the following system of equations: 

2|K]-ffl 2 [M])(A)+3[4M)M={F(rt! (40) 

Inserting equations (22-28) into equation (40) leads 
to: 

2 (41) 

Equation (41) can be written using the tensor 
notation as: 



+a,l* + ^a i a j a k b , jjkr = f- 



(42) 



This system of equations is similar to that obtained 
in equation (30) in the free case, in which i varied 
from 2 to n because the first contribution al was 
assigned. In equation (42), all of the n equations 
have a right hand side representing the forcing term 
fi*. 

In equation (42), if only one mode is assumed, and 
V * -F 

ir for 1 ^ r \s assumed to be negligible compared 

V 

with rr for simplicity one obtains: 



( 7* q ~ % 

k u -co m n )q l +-b nn a l = /, 



(43) 



* T * J * 

in which m n, 11 and 1111 are the mass rigidity, 
and non-linear rigidity terms corresponding to the 
first function respectively. 

Putting 

* 2 k\\ 

CO L =— — 

m n (44) 



Equation (43) can be written as: 



co 






J 



9 k* k* a 

^ IX j j IX i j Ci'j 



(45) 



3. One dimension non-linear frequency 
response functions 

3.1 Plate studied 

Composite laminated plates are widely used in 
military aircrafts, particularly in smaller size planes 
[30], 

The composite laminated plate (Graphite/epoxy) 
studied in reference [15] has been used in this work 
in order to compare the results obtained by the 
present model with those published in the same 
reference. The plate has five layers, with the 
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following orientation of principal axes (45°, -45°, 45°, 
-45°, 45°). Its dimensions are a = 300 mm, b= 300 
mm, h = 1 mm, and each layer has the material 
properties E L = 173 GN/m 2 , E T = El/15.4 GN/m 2 , G LT 
= 0.79 E t GN/m 2 , ht = 0.3, p = 1540 kg/m3. 

3.2 Numerical results 



The single mode assumption reduces the multi- 
degree-of-freedom system to a single degree of 
freedom one. It has been shown in previous studies 
that such an assumption may not be very accurate 
due to the increase of curvature near the clamps in 
the case of clamped-clamped beams [26], But, as 
will be shown later, the single mode approach 
remains a good tool for the estimation of the non- 
linear frequency versus the amplitude of vibration, 
in regions far from the clamps, and close to the 
centre. By assuming only one mode, the forced 
Duffing equation is obtained; the simplicity of this 
equation and the good estimation of the 
parameters of interest, lead one to use this 
approach in spite of the limitation mentioned above. 
Computing the coefficients of equation (45) for the 
composite plate for the same dimensions leads to 
the analytical formulation as: 



oo 



V 






= 1 + 3 . 267 27 9 ^ + 0 . 00694 1 46 — 



(46) 



and al is the contribution coefficient to the plate 
deflection of the first plate function. 

Equation (46) corresponding to the 1-D non-linear 
frequency response, for fully clamped laminated 
plates is represented in Figures 3-4. Figure 3 
represents a comparison between results obtained 
by the present model with those obtained using the 
hierarchical finite element method [15]. The fully 
clamped laminated composite plate is excited by a 
distributed force with a value of 4 N/m 2 amplitude. 
Under these conditions, the single mode approach 
gives a reasonable approximation for the amplitude 
of vibration of the plate at point (x*, y*) = (0.5, 0.5). 
In Figure 4, the comparison gives a bad 
approximation for the amplitude of vibration of the 
plate at point (x*, y*) = (0.75, 0.75) when the 
amplitude ratio exceeds the value of 0.1. To remedy 
to this situation, and in order to increase the 
accuracy of the frequency-amplitude relationships, 
the multimode method used for the free vibration in 
[6] has been extended to the forced case. The 
method consists of taking into account all the co- 
ordinates instead of only the single resonant co- 
ordinate and the equation obtained is a 
multidimensional Duffing equation. It can be noted 
here that the model used in reference [15] includes 
the in-plane, out off-plane functions and the terms 
containing the three first harmonics. 



In order to compare our results with those 

published in reference [15], we put ai = (Wi/h)/w* 

* 

w 

where ‘ s' are the fully clamped rectangular plate 
functions defined as 



«*•,/) =T// (*•)//(/) 

L T 

/ 

normalisation scaling factor 

G = JlAf:'(x’)f;'(y‘» 2 dx‘dy 



and G is a 
given by 

, and the 



f a 

chosen basic functions Jl are the linear clamped- 

clamped beam functions, in which V/ , for i = 1, 2, ... 
n, are the eigenvalue parameters for a clamped- 
clamped beam. 



ft (*)= 



ch(v i Yj) - cos (v. jj) sh(v i Yj) ~ sin(v. Y) 



chv. -cosy,- 



shv. -sinv. 



The values of the parameters * may be computed 

by solving numerically the transcendental equation 

coshv. cosv- =1 ... /u . ... . 

* 1 . Wi/h is the amplitude ratio at 

the centre of the plate for the fundamental mode, 



4. Multi mode approach: Numerical 

details and results 

4.1. Numerical details 

The problem consists of solving the non-linear 
algebraic system (42). This system has been solved 
numerically using the Harwell library routine NS01A. 
This routine is based on a hybrid iteration method 
combining the step descent and Newton's method, 
which do not require a very good initial estimate of 
the solution [31]. The step procedure adopted 
consists on varying the frequency parameter, which 
allowed solutions to be obtained with quite a small 
number of iterations (an average of 104 iterations 
for 18 equations). The fully clamped laminated 
composite plate is excited by a harmonic distributed 
force with an amplitude of 4 N/m 2 . The results were 
calculated using eighteen plate functions. 

4.2 Numerical results and discussion 
In order to estimate the accuracy of the results 
obtained by the present theory, and the effects of 
the approximations adopted, a comparison has 
been made with previous results, for large forced 
vibration amplitudes, based on the HFEM presented 
in [16]. 
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The co n i/ co i calculated from the present model 
appear to be very close to the results produced by 
the HFEM which used the in-plane and the three 
harmonics in the model. Figures 5-9 show a 
reasonable agreement between results given in 
reference [16] and the present results for the fully 
clamped CFRP rectangular plate. In figure 5, the 
middle-plane in plane displacements are not 
considered (pj = 0) in reference [16]; the values of 
the frequency response functions produced by the 
present model and that obtained by using HFEM are 
very similar. The comparison with the same values 
obtained by the model with that published in 
reference [16] considering the in-plane 
displacements and plotted in figure 6 show a 
slightly difference. Under the same conditions (pi = 
10, p 0 = 5), three harmonics have been used in the 
model described in reference [16], the comparison 
with the results obtained by the present model in 
figure 7 show no difference compared to figure 6. 
This has been explained in reference [16] that the 
one harmonic approximation gives a reasonable 
approximation for the amplitude of vibration of the 
plate at point (x*, y*) = (0, 0). However the interest 
has been focused near the clamps because the non- 
linearity is very important. Figure 8 shows the 
comparison between the transverse displacement 
versus co n /co*i obtained by the present model and 
that published in ref. [16] of the fully clamped 
rectangular laminated plate at the point (x*, y*) = 
(0.75, 0.75). It can be noticed here that there is no 
difference between the values. For the same plate 
and conditions, the in-plane displacements have 
been introduced (pj = 10) in the model of HFEM [16], 
the resonance frequency increases more with the 
first harmonic amplitude compared with figure 8. 
Now, when the three harmonics have been taken 
into account and introduced in the model of HFEM 
[16], the results obtained have been compared with 
that calculated using the present model and plotted 
in figure 9. It can be seen very clearly that good 
agreement has been found. 

This good agreement shows that the assumption for 
neglecting U and V in the present model can lead to 
reasonable estimation of the forced dynamic 
behaviour of composite plates. 

5. Conclusion 

A model for geometrical non-linear, steady state, 
periodic forced multimode vibration of composite 
laminated plates has been developed by using the 
Lagrange's equations, and the harmonic balance 
method. In this work, a single mode approach has 



been adopted. The results obtained indicate that 
this assumption can lead to good results in the 
centre of the plate. The accuracy of the model 
decreases near the clamps. To remedy to this 
situation, the multimode approach method has 
been used for increasing the accuracy of the 
frequency-amplitude relationships. The method 
consists of taking into account all the co-ordinates 
instead of only the single resonant co-ordinate. The 
equation obtained is a multidimensional Duffing 
equation, similar to that obtained for beams in [31], 
and for isotropic plates in [32]. The non-linear 
algebraic system has been resolved numerically and 
the results have been compared to that obtained by 
using HFEM. The comparison has been shown a 
good agreement especially near the clamps. 
Another point that should be made is that the 
model using HFEM in reference [16] take into 
account the in-plane displacements and the three 
harmonics; and the present model the in-plane 
displacements have been neglected and only one 
harmonic has been taken into account, these 
assumptions have extensively discussed in earlier 
publications [6 to 10, 17, 22 to 27]. Therefore, the 
good agreement between the results of the present 
model and that published indicates that the present 
model not only can successfully be applied to 
geometrically non-linear free analyses of 
symmetrically laminated plates, but can also 
produce accurate estimation and results of the non- 
linear forced behaviour of laminated plates. A 
further investigation considering the in-plane 
displacements and the three harmonics in the 
model is needed for more accurate results and a 
better understanding of geometrically non-linear 
dynamic behaviour. 
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Figure 3 Comparison of the forced response curve obtained by present results and that 
Published in reference [15]. Plate 2 (x*,y*) = (0. 5, 0. 5) and F d =4 N/m 2 




Figure 1: Plate notation 




Figure 2: Geometry of n layered 
laminated plate. 




<Q„1 /(O, 

Figure 4. Comparison of the forced response curve obtained 
by present results and that published in reference [15] . 
Plate , (x*, y*) = (0.75, 0.75) and F d = 4 N/m 2 
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Hgure 5. Comparison of the forced response curve obtained by multi-mode model and that 
published in reference [16] . (x , y ) = (0.5, 0.5) and F d = 4N/m 2 . 




Figure 6. Comparison of the forced response curve obtained by multi-mode model and that 
published in reference [16]. (x*, y*) = (0.5, 0.5) and F d = 4N/m 2 . 



■ multi-mode model 

O values read from graph from reference [16] 
p =5, p=10 and ussing three harmonics 

1,5 - 




Figure 7. Comparison of the forced response curve obtained by multi-mode model and that 
published in reference [16]. (x*, y*) = (0.5, 0.5) and F d = 4N/m 2 . 




published in reference [16]. (x*, y*) = (0.75, 0.75) and F d = 4N/m 2 . 
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Figure 9: Comparison of the forced response curve 
obtained by multi-mode model and that published 
in reference [16]. (x*,y ) = (0.75, 0.75) and Fd = 4 
N/m 2 
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Abstract 

In this paper, a robust adaptive control algorithm for 
speed sensorless induction motor drive is presented. 
The main feature of this control algorithm is that min- 
imum synthesis is required to implement the proposed 
strategy. It is known as model reference adaptive con- 
troller (MRAC). The MRAC appeared to be robust in 
the fact of totally unknown plant dynamics, external 
disturbances and parameter variations. The induc- 
tion motor stator resistance variations effect on the 
rotor speed are rejected by the MRAC controller due 
to its robustness. However, the performance of the 
Indirect Stator Field Oriented Control (ISFOC) are 
shading off. It is so because the reference commands 
computation of the ISFOC are based on the induction 
motor model, which includes that resistance. There- 
fore, a combination of three observers is used : a Lu- 
enberger observer to obtain the rotor flux and a two 
MRAS observers for rotor speed and stator resistance 
estimation. Digital simulation results are presented 
to highlight the improvement in performances of the 
overall proposed control scheme. 

Keywords: Induction motor, ISFOC, MRAC, 
sensorless rotor speed, stator resistance, MRAS and 
Luenberger observer. 



1 Introduction 

In concert with the advances in factory automation, 
servo systems have become indispensable to many ap- 
plications in industry. Process in the use of induction 
motors with high performance AC servo systems is 
remarkable. Previsiouly, DC motors were used ex- 
tensively for variable speed drive because of ease of 
control. However, AC motors have been controlled 
like DC motors using mainly a field oriented control 
approach [1]. Because of advances in power electron- 
ics and microprocessors, the induction motor drive 



used in variable speed control has become more at- 
tractive. In view of the advances in power electron- 
ics and microprocessors, digitally controlled induction 
motor drives become increasingly popular. In many 
industrial drives, advanced digital control strategies 
for the field oriented control of induction motor with a 
conventional speed controller like such as the PI con- 
troller, have gained wide acceptance in high perfor- 
mance AC servo systems. 

A high performance servo system must have good dy- 
namic speed command tracking and load regulating 
responses, and the performance must be insensitive 
to system parameter variations. The conventional PI 
controller is easily implemented and highly effective if 
the load changes are small and the operating condi- 
tions do not take the system too far from the linear 
equilibrium point. However, in many applications as 
machine tools, the drive operates under a wide range 
of load changes and the system parameters vary sub- 
stantially. In addition, because, it processes only one 
degree of freedom, the PI controller can’t provide good 
tracking and regulating performance simultaneously. 
To overcome this drawback, the control algorithm 
should include a complicated computation process to 
eliminate the variations in load disturbance and sys- 
tem parameters and also obtain high performance AC 
servo systems. However, the control algorithms which 
are suitable to these systems have become increasingly 
more complicated, requiring extensive computations 
for real time implementation [2]. 

Hqq loop shaping is used in the speed loop [3] and 
[4]. It remains robust against disturbances and load 
torque. However, this approach requires the knowl- 
edge of the limits of the disturbances. The sliding 
mode mode control is also used in the speed loop [5] 
but the most drawback of this type of controller is its 
higher switching frequency [6]. In recent years, speed 
controllers based on artificial intelligence techniques 
such as fuzzy logic [2], [7], [8] and [9] and neural net- 
work have [10] been widely proposed. 
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Since, these approaches don’t require the knowledge 
of a mathematical machine model, the algorithm re- 
mains robust despite parameter deviations and noise 
measurement. However, the computation expenses 
and the requirement of expert knowledge for the sys- 
tem set up have seriously restricted their applications 
in practice [2]. 

Parameter variations represent a great problem not 
only for the speed controller, but also in the ISFOC 
decoupling of the induction motor [11]. As, the decou- 
pling equations are deduced from the induction motor 
model which depends on the induction motor param- 
eters. To overcome this problem, the identification 
and on line tuning of these parameters are strongly 
needed. 

All high performance vector controlled induction mo- 
tor drives depend on the quality of the signal measure- 
ment specially the rotor speed. As, it is used either in 
the speed control loop as in the field weakening oper- 
ating region [12]. The principle of the field weakening 
operation consists of varying the flux reference in pro- 
portion to the reverse of the speed. 

In a published research work [13], a direct adaptive 
control law is used in the speed loop. However, the 
effectiveness of the proposed method is restricted only 
to the load torque. Hence, robustness against motor 
parameter Variations as a major factor are not con- 
sidered. 

In this paper, a direct adaptive control based on the 
Model Reference Adaptive Control (MRAC) is per- 
formed for the speed control loop. The effectiveness 
of this algorithm is tested both with load torque and 
stator resistance variations. The ISFOC is used to 
assume the decoupling of the induction motor, the 
variation of stator resistance makes fails for this type 
of decoupling. To compensate the effect of the stator 
resistance variation, an algorithm based on the Model 
Reference Adaptive System MRAS principle is pro- 
posed to make an online identification and decoupling 
process tuning. The rotor speed is reconstructed by an 
MRAS speed observer. Digital simulation results are 
provided to show the improvement in performance of 
the proposed robust control of a sensorless IM drives. 
The paper discusses the Adaptive Control for Sensor- 
less Induction Motor Drives. Section 2 presents stator 
flux orientation model. Section 3 illustrates Speed and 
stator resistance adaptive rotor flux observers. Sec- 
tion 4 introduces the adaptive scheme for stator resis- 
tance and rotor speed observation. Section 5 shows 
the design of an adaptive speed controller. Section 
6 explains the simulation results and section 7 is the 
conclusions. 



2 Stator flux orientation model 



4*ds — 0 s and (j)qs — 0. 

Assuming these conditions, the dynamic model of an 
induction motor can be represented according to the 
d and q axis components in a synchronous rotating 
frame as [11], [16], [17] and [18] : 



Vds — T's ids T $ *&ds 

Vqs T*s iqs T Ld s *&ds 
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It can be seen if the stator flux is kept constant, the 
electromagnetic torque can be only controlled by the 
q axis current. 

3 Speed and stator resistance 
adaptive rotor flux observers 

Assuming that linear magnetic circuits, equal mutual 
inductances and neglecting iron losses, if the stator 
currents and the rotor flux are selected as the state 
variables, the state equations can be expressed in a 
stationnary reference frame as [18] and [19]: 
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For ISFOC [1], [14] and [15], the stator flux vector is 
aligned with the d axis and setting the stator flux to 
be constant and equal to the rated flux, which means and 
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i s = [ ias ifis ] T is the stator currents vector, 
v s = [ v as vp s ] T is the stator voltage vector, 

4>r = [ 0ar <t>(3r ] T is the rotor flux vector. 

r s and r r are the stator and rotor resistances 
respectively, L s , L r and M are the stator and rotor 
self and mutual inductances respectively, a is the 
leakage inductance coifficient and cu is the rotor 
speed. 

The electromagnetic torque developed by the motor 
is expressed in terms of rotor flux and stator currents 
as: 

T e =p* $ r (gi i s , (8) 

While the load torque acts as a disturbance ac- 
cording to the mechanical equation: 

^=T,-T L . (9) 

mi is the rotor inertia, p the number of pole pairs and 
Tl the load torque. 

The adaptive rotor flux observer, to get the estimated 
stator current and rotor flux using mainly the mea- 
sured stator current i s and the stator voltage v s , is 
described by the following state of equations: 



an adaptive mechanism to turn out the observed vari- 
able. 

In this work, the actual system is considered as the 
model reference and the observer is used as an the ad- 
justable one. 

The stator resistance [22] and [23] and the rotor speed 
[24] are reconstructed by using the following adaptive 
mechanisms: 

t 

v se = Kp\c^ s i se T Kji J c^ s i se dt (H) 
0 
t 

LU e — KpCis^re T K j J & is(firedt (12) 

0 

Where Ki , Kn, K p and K p \ are the integral and pro- 
portional constants respectively. As shown in figure 
1 , 7Q, Aii, Kp and K p \ are considered as the weights 
multiplying the product between current errors and 
the stator currents and the estimated rotor flux and 
its integral respectively. Therefore, the tracking per- 
formance of the speed estimation and the sensivity to 
noise are depending on the proportional and integral 
coifficient gains. The integral gain Ki is chosen to be 
high for fast tracking of speed. While, a low propor- 
tional Kp gain is needed to attenuate high frequency 
signals denoted as noises [17]. 

MRAS representation structure for identifying the ro- 
tor speed and the stator resistance simultaneously is 
given in figure 1. 



A 


i se 




Tllle 


7ll2e 




^ se 




' B 1 ' 


dt 


(j) r e 




A21 


^22e 




0re 


+ 


. B 2 . 



Le(t ) 



(pre 

( 10 ) 

Where , ”e” , denotes the estimated values, L is the 
feedback gain matrix of the Luenberger observer and 
e(t) is the stator current observation errors. 

In this work, the matrix gain observer is obtained us- 
ing the principle of the LRQ as proposed in [15], [16] 
and [18]. 



4 Adaptive scheme for stator re- 
sistance and rotor speed ob- 
servation 

The rotor speed and the stator resistance are recon- 
structed using the model reference adaptive system 
(MRAS) [18], [20] and [21]. The MRAS principle is 
based on the comparaison of the outputs of two es- 
timators. The first is independent of the observed 
variable named as model reference. The second is the 
ajustable one. The error between the two models feeds 




Figure 1: MRAS observed rotor speed and stator re- 
sistance. 



5 Design of adaptive speed con- 
troller 

The dependance of the induction motor parameters 
on the temperature is a serious hindrance in high per- 
formance control of the induction motor drives [5] and 
[25] . The temperature rise is caused by power losses 
in the motor [26] specially when the stator current ex- 
ceeds it nameplate level. Under certain circumstances, 
such as overload or start up operations, the stator cur- 
rent can be much higher than the rated one, resulting 
in a significant increase of the stator resistance [27]. 
Compared to a cold motor, the induction motor pa- 
rameters can increase by an order to 100% [25] and 
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[28]. The use of the classical controllers such as the 
proportionnel and integral controller is insufficient to 
provide good speed tracking performance. To over- 
come these problems, an adaptive controller based on 
the MRAC principle is proposed for the speed control 
law synthesis [13] and [29]. 



After replacing u(k) with its expression in (13), the 
model output can be written in the following form: 



Vp{k) 



bjt 0 (k)q 1 

1 + («1 + biSo(k))q~ 1 



r(k ) 



(20) 



5.1 Discrete model reference adaptive 
control 



we consider the single input single output discret lin- 
ear time invariant plant of the first order system. It 
is described by the input - output pair {u(k), y p (k)} 
and defined with the transfer function G(q~ 1 ) by: 



rln -u = y_pW = M 1 

yq 1 u(k ) 1 + ai< 7 -1 



(13) 



A model reference represents the desired behavior of 
the system. This model is described by the input - 
output pair {r(fc), y m (k)}. It is identified by the fol- 
lowing transfer function: 



Gm{q x ) 



Vm (k) 

r(k) 



bmQ 

1 tflmf 1 



(14) 



Where: 

r(k) is the reference input, 

ym{k) is the output of the model reference. 

The reference input r(k) is bounded. 

The objective of the numerical model reference con- 
trol law u(k) is to minimize the tracking error. This 
control law satisfies the following equation [30]: 



lim e(k) = lim (y p (k) - y m (k)) = 0 (15) 

k — >-oo k — >-oo 

The minimization criterion is defined by: 

J = \e(k) = ± (y p (k) - y m {k )) (16) 

We designate by 0 the vector of the controller param- 
eters. To update these parameters, we adopte the 
gradient algorithm. It is defined by: 

e{k + i) = e{k)-a^ ) (i7) 

With: 

a is the adaptation factor. 

The control law which minimizes the criterion J is 
given by [31]: 



u{k) = t 0 (k)r(k ) - s 0 (k)y p (k) (18) 



Where t Q (k) and s Q {k) are the components of the con- 
troller vector. 

9{k)= [ t Q (k) s Q (k ) ] T 

By replacing each element with its expression in (17), 
we obtain: 



f + 1) 

\ s o{k + 1 ) 



o(k) \ 

o(k) J 



de(k ) 
dt 0 (k) 
de(k ) 
ds 0 (k) 



(19) 



To formulate the error, we subtract the reference out- 
put y p (k) of (20) to model output. Subtracting y p (k) 
from y m (k) of (14), we obtain: 



e(k) = 



b\t 0 {k)q 1 



b m q 



-1 



1 + (ai + bis Q (k)) q 1 1 + a m q 



-l 



r(k ) 

(21) 

The differentiation of the error through 0{k) can be 
written in the following form: 



de(k) 

d6(k) 



l-\-(ai-\-b 1 s 0 ^k))q 1 r (k) 



(22) 



In order to obtain a dynamic system equivalent to the 
reference model, we impose the characteristic equa- 
tion. Therefore, 



de{k) 

d6(k) 



big 1 

1 ~\~ a mQ 1 
big -1 



r(k) 

Vp{k) 



(23) 



Repiacing ^ 



with its expression in (19), we obtain: 



I to (k + l)=t 0 (k)-a J ^ T r(k) 

[ s 0 (k+ 1) = s„(k) + a 1 ^ g _ 1 y p (k) 

The developpement of these equations yields: 

to(k) = (l-a m )t 0 (k-l) +a m t 0 (k-2)-ae(k)b 1 r(k) 

(25) 

so(k) — (1 dm) so(k 1) 4~ so(k 2)Ttr e(A:) b\ y p (k^ 

(26) 



The structure of the direct adaptive control with the 
MRAC principle is shown in figure 2. 




Figure 2: The basic structure of the Model Reference 
Adaptive Control. 
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5.2 MR AC for induction motor drive 



If we suppose that the load torque is as a disturbance, 
then using the mechanical equation, we obtain: 



c j(s) k 

T e (s) l + TS 



(27) 



where k = — and t = — 

mi mi 

The discrete transfert function can be written as: 



uj(k) bi q 1 
T e (k) 1 + cq q- 1 



(28) 



The model reference must be chosen for the same 
plant as the real one. Therefore, 



Um{k) _ b rn q 1 
^ref{k') I T Q 



(29) 



Where b m and a m are the model reference parameters. 
The command law of the model reference is given by 
the recurrent equations: 

u(k) = t 0 (k) u ref (k) - s 0 (k) u e (k) (30) 

to(k) and so(fc) are updated from the two following 
recurrent equations: 

t 0 (k) = (1 - a m ) t 0 (k - 1) + a m t 0 (k - 2)-ae(k) 
bicj ref (k - 1) 

(31) 

s 0 (k) = (1 - a m ) s 0 (k - 1) + a m s 0 (k - 2)+cm(k) 

blCJref (k - 1) 

(32) 

The proposed configuration of the robust adaptive 
control law is detailed in figure 3. 




Figure 3: Configuration of the robust adaptive sen- 
sorless control system. 



6 Simulation results 

The dynamic behavior of the system was investigated 
using computer simulations with Matlab / Simulink 



software. 

We have tested the robust adaptive control algorithm 
for sensorless induction motor at different speed tar- 
gets under no load and with load applied. The effects 
of stator resistance variations is also studied. 

The results obtained at variable step speed targets un- 
der no load are reported in figures 4 to 7. 

It can be clearly seen that the d and q axis stator 
flux converge to the nominal value and to zero re- 
spectively. This indicates that the decoupling of the 
induction motor is established and ensured. 

The system and the model reference responses at no 
load operation are given in figure 6. It can be seen 
that a good agreement between system and model 
reference responses is obtained. The feedforward and 
feedback gains are shown in figures 4 and 5 respec- 
tively. It is noticed that when an error between the 
system response and the model reference one occurs, 
the feedforward and feedback gains react to minimize 
this error. They converge to a fixed values when the 
error between the system and the model reference re- 
sponses becomes zero. 

To investigate the behavior of the proposed sensorless 
control drives against external disturbances, a load 
torque was added. The results obtained under load 
applied are shown in figures 8 and 9. The obtained re- 
sults show that an error between the system response 
and model reference one occurs. It converges to zero 
thanks to the on line adaptation scheme of the con- 
troller parameters. 

The robustness of the proposed sensorless control 
drives is tested according to the stator resistance vari- 
ations. The obtained results are given in figures 10 
and 11. It is shown that the system response converge 
to the model reference response. An error between 
the system and the model responses occurs when sta- 
tor resistance varies. It is compensated thanks to the 
adaptive identification of the controller parameters. 
As it is seen in figure 11, a decoupling error occurs at 
a time of the stator resistance variation. But, this er- 
ror is maintained and there is no compensation. This 
is because of the ISFOC commande equations which 
are established from the nominal value of the stator 
resistance. 

To overcome the sensivity of the vector control algo- 
rithm (ISFOC) against stator resistance variations, an 
identification and on line adjustment of the stator re- 
sistance in ISFOC algorithm is necessary to maintain 
the decoupling of the induction motor. We have used 
the MRAS principle to observe the stator resistance 
value. The obtained results with on line adjustment 
of the stator resistance are reported in figure 13 to 15. 
The system response converge to the model reference 
one. An error occurs at a time of the stator resistance 
variation. It converges to zero thanks to the adaptive 
control algorithm. The target, observed and actual 
components of the stator flux are shown in figure 14. 
The d and q axis stator flux converge to the nominal 
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value and to zero respectively despite of the stator 
resistance variations. This is because of the on line 
identification and tuning of stator resistance in the 
ISFOC algorithm. 

The observed and the actual stator resistance are 
given in figure 15. The results show that the observed 
stator resistance follow exactly the actual stator re- 
sistance, indicating the validity of the proposed algo- 
rithm. 




Figure 4: Evolution of to(k) at variable target rotor 
speeds under no load. 





Figure 7: Actual, observed and target d - q compo- 
nents of stator flux at variable target rotor speeds un- 
der no load. 
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Figure 8: Actual, MRAS observed and output of the 
reference model at variable target rotor speeds under 
load. 



Figure 5: Evolution of so(k) at variable target rotor 
speeds under no load. 
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Figure 6: Actual, MRAS observed and output of the 
reference model at variable target rotor speeds under 
no load. 




Figure 9: Real, observed and target d - q components 
of stator flux at variable target rotor speeds under 
load. 
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Figure 10: Actual, MR AS observed and target rotor 
speeds without on line adjustment of stator resistance. 
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Figure 13: Actual, MRAS observed and target rotor 
speeds with on line adjustment of stator resistance. 
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Figure 11: Actual, observed and target d - q com- 
ponents of stator flux without on line adjustment of 

stator resistance. 
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Figure 14: Actual, observed and target d - q compo- 
nents of stator flux with on line adjustment of stator 

resistance. 
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Figure 12: Actual stator resistance trajectory varia- 
tion. 



Figure 15: Actual and MRAS observed stator resis- 
tance. 
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7 Conclusion 

In this paper, we have presented the synthesis of a ro- 
bust adaptive control algorithm for sensorless induc- 
tion motor drives . The proposed method has been 
applied to an indirect stator field oriented induction 
motor drive system. The adaptive with control al- 
gorithm for the rotor speed shows good performance 
under no load, load applied and against stator resis- 
tance variations. The influence of the stator resistance 
variation on the ISFOC algorithm was eliminated by 
the proposed algorithm for the identification and on 
line ajustement of the stator resistance. The validity 
of the robust adaptive control algorithm for sensorless 
induction motor drives has been verified by simulation 
under Matlab/Simulink environment. 

Practical implementation of the proposed method is 
a subject of future follow up research work. 
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Abstract 

This paper will provide an effective help to all those 
who make decisions regarding the planning and 
implementation of projects wind Energies. We used 
the software Alwin for sizing a wind farm of 140 MW 
installed in Melloussa cite, 25 Km from Tangier, in 
north Morocco. We analyzed the collected wind 
data during one year, to assess wind power 
potential, the prediction of produced electric power 
in the site based on a judicious choice of Wind 
Turbine Generator, as well as, the automatic 
determination of direction of the wind on the site 
and analysis of wind turbulence. Thus, the 
integration into the grid provides several benefits. 
However, this type of production which has an 
impact on the electric network and its connection is 
subject to technical constraints. We will study the 
various generators used to produce electric power 
from wind energy and choose the one with a good 
price per kilowatt hour, then to study the feasibility 
of integration of this wind farm in the electric grid. 
And finally, we treated the various principles, 
concepts used in data transmission and remote 
monitoring of wind farms, and we proposed 
industrial solutions of telemanagement which can 
be applied in remote wind farm in the site. 

Keywords: Wind data, Wind park, ALWIN software, 
Aerogenerators, Flicker, harmonics, Telemanagement. 



1. Introduction and the general 
objective 

To optimize the exploitation of the wind farms, 
some measures must be taken into consideration. 
Indeed, a bad choice of certain parameters could 
harm a wind turbine establishment. The reason way 
a pre-study phase appears of extreme importance 
to identify the conditions and constraints to be 
taken into account in the realization of a wind 
project [l]-[2]. 

We must therefore proceed with a study and 
analysis of the climatology of the site considered, 
since a good knowledge of the site allows a better 
exploitation and production of electrical energy. 
However, the coupling of the Wind turbines at 
national grid has an impact on power quality and 
dependability. As a result, the connection is subject 
to certain technics and requires some adjustments 
in the network to ensure a proper functioning of this 
latter. 

Thus, to better optimize and secure wind farms, a 
control system and remote monitoring must be 
established. This latter must be designed to allow 
automatic operation of wind farms. 

And therefore, the general objective of this work is 
to help develop indigenous energy sources to 
supply electricity and profit from great potentials of 
wind power in Morocco as long it’s a clean, safe and 
a free source of energy. 

This paper is organized as follows: section 2 explains the 
Study feasibility of dimension of wind farm via Alwin 
software, and at section 3, we present the results of 
dimensioning a wind farm of Tangier, in section 4 we 
introduce Impact of the production of wind farm on the 
national network, we describe in section 5 the Concept of 
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remote wind farm. In the end, section 6 provides 
conclusions and directions for future research. 

2. Dimension of wind farm via software 
Alwin 

2.1. Sizing techniques: [3] -[4] -[5] -[6] 

2. 1.1. Wind power: 



p = ( y 2 )psv 3 


(1) 


Vmoy = Y, Vift 


(2) 



P: wind power (kw), 

P the density of air (Kg/m 3 ) 

V: the average wind speed in (m / s) 

S: the area swept by the wind turbine in ( m 2 ) 

Vmoy: the average wind speed in (m / s) 

fi: Frequency in% of the wind speed Vi in 

(m/s)to class i. 
i : class of wind speed 

Vi: the wind speed in (m / s) at fi frequency 

2.1.2. Energy provided by a wind turbine! 

E = T t x’Z p Mf l <y) ,,, 



With: 



E: 


energy produced kWh. 


Th: 


time in hours 


Pi (V): 


power (KW) given by the wind speed Vi 


fi(V): 


frequency of wind speed Vi 



2.1.3. Capacity Factor 

The capacity factor assesses the degree of utilization 



N2: number of row of turbines 

N: total number of turbine to be placed on the 

site 



Mathematical modelling of frequency distribution of 
wind: distribution WEIBULL: 

As it is difficult to manipulate all of data on a 
frequency distribution of wind, it is more convenient 
for theoretical considerations to model the 
frequency histogram of wind speeds by a 
mathematical function that continues through 
discrete values table. We can therefore choose the 
model WEIBULL [3]-[4]-[5]. Indeed, for periods 
ranging from weeks to a year, the WEIBULL function 
reasonably represents the observed rates. There is a 
probability density function in the following form: 

P(V) = (K/C)(V/Cf- ] expHV/c/) (6) 

With: 

P (V): probability density of the velocity V 

K form factor of the curve (dimensionless) 

C: scale factor of the curve (meter / second) 

The average wind speed can be found by 
integrating the probability density is therefore the 
formula (6): 

V = \VP(V)dv (7) 

moy J w 

Or =CT(\ + l/k) (8) 



Where T (x) is the mathematical function Gamma, 
defined by: 



/•CO 1 

r(jc) = f e^t^dt 

JO 



(9) 



of the turbine, it is defined as follows! 

«%M100 <P„)/P m 



(4) 



These two powers Pmoy and Pmax are chosen for 
the wind turbine. As the capacity factor is high, it is 
certainly good choice of the wind turbine. 

Pmax: maximum power delivered by 

aerogenerator 

Pmoy: average power delivered by 

aerogenerator 

2.1.4. The total number of turbines to be placed on 
the site [1] Conditions to be attained: 

(N1 +1)* 10H <1 
(N2 +1)* 3D <L. 

N = N2 * N1 (5) 

Key: 

I: dimension of the field perpendicular to the 

prevailing wind direction. 

L: size of land parallel to the prevailing wind 

direction 

D: rotor diameter of the machine 

H: height of the tower 

Nl: number of turbines per row 



We can deduce other important parameters such as 
the average power: 

P^=0.5pC 3 T(l + 3/k) (10) 

Thus, the distribution WEIBULL can facilitate many 
calculations that are necessary for the analysis of 
wind data. 

2.2.1. Determination of WEIBULL parameters [4]-[5]: 
Now, we must focus our attention on the 
parameters K and C, as they are the key to using the 
WEIBULL distribution. These parameters define the 
shape and scale of the curve WEIBULL. 

The form factor, K (dimensionless), expressed in the 
form of WEIBULL curve: A high value of K implies a 
narrow distribution, with winds concentrating 
around a value, while a low value of K leads to 
widely winds dispersed. 

The factor C (m/s), determines the position of the 
curve. Generally, its value is high for windy sites and 
low for quiet sites. Various methods can be applied 
to determine K and C from a given wind distribution. 
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c =1.125 v mo> , /(I -B) 

K = 1 + 0.483 (V moy - 2) 051 

B =1-0.81 (V moy -l) 0 089 (n: 

This probability density function is called Rayleigh 

distribution when K = 2 

2.2.2. Changes in wind direction: 

We can determine the dominant wind direction. The 
alignment of the turbine is intended to capture the 
maximum energy available in the wind. Thus, in so 
far as the ground permitted, wind turbines should 
be aligned perpendicular to the predominant 
directions of the wind. In general, the wind direction 
is presented as a wind rose. 

The parameters in the paragraph above "sizing 
techniques" describe the characteristics of wind in 
any site. 

Winds are generally characterized by: 

• The average speeds (hourly, daily, monthly 
and annually). 

• The predominant directions. 

• The frequency of each of the speeds and 
directions. 

• Diurnal Pattern / Duration of calms / calm 
analysis. 

• A good knowledge of these characteristics 
is important in so far as they are involved 
in: 

o The choice of wind turbines and the 
evaluation of their performance; 
o The implantation (Proper orientation of 
wind). 



3. Presentation of results of 

dimensioning [5]- [7] 

3.1. Wind Potential of Tangier’s site 

We processed and analyzed the data of station in 
Tangier of 369 days for an assessment of the wind 
potential using the menu under "Annual Review". 
We found the frequency distribution of wind speed 
as follows: 




Figure 1. a: Statistic of Wind potential of Tangier’s site 



Site: Altitude: 140 m 

TANGER Temperature: 22.8 *C 

Start Meas : 01.12.00 00:00 No. of Averages: 3B9 darys (560538 d 1 min.) 



Frequency, % Wind Speed Statistics 




Figure l.b: Wind potential of Tangier’s site 



The annual review of the frequency distribution of 
the wind speed gives the following information: 

The average speed of the site is: 6.04m/s. The 
maximum speed of the site is: 25.4m/s. Speeds 
predominant site varies between 2.5m/s and 11 m/s. 



3.2. Daily variation / calm analysis [5]-[7] : 



Site: Altitude: 140 m 

TANGER Temperature: 15.0 “C 

Start Meas.: 01.03.01 00:00 No. of Averages: 31 days (44640 a 1 min.) 




Figure 2. Daily variation / curve calm. 



On the one Hand, the analysis of wind turbulence in 
the site through the curve of the daily model allows 
us if there is a continuity of service (Permanent 
Power) at the station. On the other hand, the curve 
calm analysis of the station in the form of sticks 
shown in Figure 2 shows the cumulative frequency 
percentage for speeds below 2.5 m/s and 4.5 m/s 
and remain for some time calm. The determination 
of the frequency of calm is needed to determine the 
speed boot of aerogenerator which will be used for 
energy production. Also, we need the frequency of 
calm for determining the storage capacity of 
batteries included in the electrification systems 
(independent), as well to ensure continuity servicing 
production of energy. 

3.3. Predicted annual energy produced by the wind 
farm based on the appropriate choice of wind turbine 
[5]- [7] 

The submenu "Prediction of Electric Power" is 
designed to show the user the power measured in 
MWh, the average power of the wind turbine 
selected and the capacity factor, which assesses the 
degree of utilization of wind turbine, 
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TANGER NORDEX N54/1000 

Time: 01.12.00 00:00 Meas. Height: 1 0.0 m Hub Height: 54.0 m 

Altitude: 140 m Temperature: 22.8 °C Conversion factor (1 /7-Power Law): 1.272 




Figure 3. a: Curve prediction of energy. 

Wind Energy Converter Catalog: DEWI9503.CAT 

NORDEX N54/1 000 Hub Height: 54.0 m 




Figure 3. b: the machine power curve characteristic 

- The choice "NORDEX N54/ 1000," Pn = 1MW can 
produce an energy of 3.69 GWh annually by each 
turbine. 

Knowing that the installed capacity of 140MW, or 
140 Wind turbines give the annual energy 
production of 516.6 GWh in the site. 

3.4. Result of wind Direction [5] -[7] 




Figure 4.a:Wind direction Statistic: Tangier’s site. 



SlK Aitudc: WDm 

TANGER rampin'.*? 326"C 

STSKMfa* ni. 12.00 00.00 NO. mil.) 




Figure 4.b: Wind rose Tangier’s site 
The annual review can opt for the north-east, 
specifically an inclination with respect to North of a 



varying geometric angle between the first and the 
third a sector (between 15 degrees and 45 degrees). 

4. The Impact of the wind production 

park on the national network 

4.1. Phenomena of flicker and harmonics [8]-[9] 

The flicker [8] 

The flicker is a visual sensation produced by rapid 
changes of an effective tension. 

The main cause of flicker emission by the wind 
turbines is a turbulent wind. 

- Other phenomena can cause flicker as: 

*The shadow effect of the tower. 

* The blades error pitch. 

* The rotor orientation error. 

The Harmonics 

The harmonics are introduced by elements of 
electronics power, especially the frequency 
converters, the thyristors controllers and capacitors. 

4.2. The choice of most appropriate turbines 

4.2.1. The different technologies of wind turbines 
(aerogenerator) 

Table 1 summarizes the different technologies of 
wind turbines, we opted for the DFTG for the 

following reasons: 

• Possibility to regulate the tension in 
connection point where it is injected (very 
small flicker) [ 10 ] 

• Asynchronous double fed machine gathers 
the advantages of asynchronous and 
synchronous machines. 

• Rugged and allows the inspection of the 
reactive energy. 

• Costly compared to the asynchronous 
machine, but it is less than the synchronous. 

A kilowatt hour price much better the than 
synchronous machine under the same regime of wind 
[ 10 ]. 

4.2.2. The choice of the machine double fed 
induction asynchronous generator (DF1G) 




Figure 5. The double-fed asynchronous machine 
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Table 1. Summary of the different technologies [6]-[8 ]-[9 ] of wind 
turbines 



Speed 

control 


Type of 
generator 


Reactive 

power 


Disturbance 
created in 
the network 


remark 


Fixed 

speed 

With 

gearbox 


Synchronous 


Can be 
controlled 
by the 
current 
excitation 


Mechanical 
oscillation 
because the 
generator is 
connected 
to the main 
frequency 
(flicker). 


More used in 
large units 


Asynchronous 


Requires a 
correction 
device 
cos(cp) 


Requires a 
soft start 
device. 
Flicker 
generation 


Use in the past 
in small units. 


A speed 

adaptable 

within 

limits 

with 

gearbox 


Asynchronous 
adaptation of 
rotor slipping 
resistance 


Requires a 
correction 
device 
cos(cp) 


-Soft start 

device 

-Possible 

Harmonic 

by 

adjusting 

the 

slipping. 


Poor 

performance. 


Variable 

Speed 

Rotor 

With 

gearbox 


Asynchronous 
dual power 


Power 
factor 
control by 
frequency 
converter 
rotor. 


Flicker 
reduced but 

harmonics 

Problem 


Solution 
technically 
interesting but 
we must see 
the investment 
cost 


Asynchronous 
dual rotor 
winding 


Requires a 
correction 
device 
cos(cp) 


Requires a 
system of 
soft start 
device 


Solution 
technically 
interesting but 
we must see 
the cost of 
investment 


Rotor 

without 

Variable 

speed 

gearbox 


Multiple 

Synchronous 


Can be 
controlled 
by the 
current 
excitation 


Little 

disruption 

on 

networks. 


Solution 
technically 
interesting but 
greater cost 
compared to 
the DFIG. 



4.3. The fliker emission calculation 

4.3.1. In continuous operation [9] 

p* =Pn = N 05 .C('¥).-^- (9) 

C (\| /): coefficient flicker of 10.4 (value given by the 
standard 61400-21). 

Sn is the rated apparent power of the wind turbine. 
See is the apparent power of short-circuit at point of 
common coupling (PCC). 



According to international standard: IEC 

(International electrotechnical committee 61400-21 
of measurement and evaluation of power quality 
characteristics of grid connected wind turbines. IEC 
61400-21 involves: N10 = 10 and N 120 = 120. 
Pst: Compatibility levels at short time for flicker 
networks Low Voltage and Medium Voltage. 

Pit Compatibility levels at long time for flicker 
networks Low Voltage and Medium Voltage. 

To evaluate the worst flicker, we take Kf (\j/) = 1. 

4.3.3. Presenting the results of flicker calculations at 
Pcc [9]-[ll] 



Table 2. The results of the flicker calculations at Pcc 



Active 

power 

(KW) 


Power 

(KVA) 


Numb 
er of 
unity 


Post 

connectin 

g 


Sccmin 
(MV A) 


Emission of flicker 












continious 


During 

switching 

operations 












Pst=Plt 


Pst 


Pit 


1000 


lin.ii 


140 


Wind Park 


2779 


0.0491 


0.0608 


0.0584 


Post of 
Melloussa 


3379 


0.0404 


0.05 


0.0480 


Post of 
Tanger 


3199 


0.0427 


0.0528 


0.0507 



The flicker factors below 1, which complies with the 
standard, have been obtained. 

4.4. Harmonics Calculation at PCC: 

The harmonic current order h at the PCC due to N 
wind turbine [9]: 



Key: 

Ihi: the harmonic distortion of current order h due 

to a single wind turbine. 

n: the report of a transformer wind turbine. 

The coefficient a according to the order of 
harmonics. 

The rates of harmonic tension due to 140 wind 
turbines at bus bar MT are in table 3. 



4.3.2. During operations commutation 
P sl =^[N.N w (Kf^).S n ) 3 - 2 f 31 

Pu = -N no (kf ((/)S„) 3 ' 2 ]° 31 



Kf(l//) Lactor: is the flicker of the wind for the phase 
angle impedance given network 



Table 3. The rates of harmonic tension due to 140 wind turbines 





harmonic 

pair 


odd 

harmonics 


total 

harmonic 


Value relative 
to 22kV 


0 ,28% 


1,19% 


1,22% 



The ONE (National Office of Electricity) requires the 
following rates of harmonics [ 11 ]: 
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kind of harmonic 


Pair 


odd 


Rate of harmonic 


0,6% 


1,5% 



• The rate of overall harmonic is therefore 
well below the 1.6% imposed by ONE [11]. 

• The rate of harmonic issued by all the wind 
turbines is not embarrassing for the 
network and does not exceed the norms 
(limits prescriptives). 

5. Concept of remote wind park 

5.1. Purpose of the remote [12] 

The remote is a powerful tool in management and 
control. It improves the speed of intervention of 
maintenance services and increases the 
performance and availability of machines. It is 
transmitted on a physical data and messages of 
breakdowns to a room monitoring. 

5.2. Solutions proposed for remote for the park of 
Tangier 

5.2.1. Problem Analysis 

We want to monitor the wind farm in two different 
locations: the village post of Melloussa and the post 
of the city of Tangier as shown in the following 
diagram: 



900m 25 KM 




5.2.2 Data Transmission 
Possible solutions 

Each compact station at the foot wind turbine has a 
telephone line from the central Melloussa. The latter 
is linked to the central Tangier by a high voltage 
(225kV), we adopted the following solution: 




- RTC (Reseau Telephone Comute: phone network 
commuter ) and transmission power line on the high 
voltage line 

Forwarded by RTC and current holders on the high-voltage line: 

- Current holder link 




I/E : injector/extractor. 



- Reasons for the choice: 

It’s reliable, economic (cable transmission is the 
same cable HV) and already operated by ONE 
(National Office of Electricity) in Morocco at the 
factory IDRISS 1. 

5.3. Control and monitoring of wind turbines [5] 

5.3.1. Solution based on PLC (Processing logic 
control) 

The configuration of the wind farm permits to 
consider only one solution if it is based controller, 
for each turbine, we will install a mini-controller. 

5.3.2. Reasons for choice: 

Maximum automation for a minimal investment, 
using a stand-alone and far simplicity of installation, 
programming and manipulation [13]-[14]-[15]. 

5.3.3. Architecture of the solution to be applied in 
the park of Tangier 

At the foot of each turbine mini-PLC is installed. 

The CPU will be installed in the room Melloussa and 
manages all micro-controllers before contacting the 
supervisor PC installed in Tanger. 




Figure 6. Architecture of the solution to be applied in the park of 
Tangier 

CPU: Central programmable unity. 

PLC: Processing logic control (PLC micro) 

The micro-robot (SIMATIC S7-200) meets the 
desired functionality [9]. Already operated by the 
wind Park of Tetouan in Morocco, the only 
drawback its high cost. 

5.3.4. Solution based on a microcontroller 

The number of entry / exit and the size of the 
programs at each turbine are not important, so a 
microcontroller based solution seems very 
reasonable, because it replaces the role of the 
controller with a lower price. 

5.3.5. Choice of network topology [13]-[14]-[15]: 

In the star topology, all stations of wind turbines are 
connected directly to a server or a hub that is the 
central node through which all transmissions pass. 
This topology allows to easily add equipment (one 
cable per device in the limit of the server capacity). 
Network management is facilitated by the fact that 
the equipments are directly searchable by the server 
and all transmissions spend (centralization of 
software). Moreover, a failure of terminal equipment 
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does not affect the operation of the rest of the 
network. However, it can lead to significant lengths 
of cable. It is more economical than the parallel bus. 
And given the network structure, this solution is very 
appropriate for our application. 

In the case of topology bus, if the device is the 
optical type, the adding operation is more difficult 
because it requires cutting the fiber to the location 
of the connection and that the stations located 
downstream report to the server will be blocked. 
NOTE. - The connection between the PC Melloussa 
and wind turbines may also be performed by a RS 
485: 

There are two types of interface circuits: 
asymmetrical and differential. 

5.3.6. The asymmetric-RS232 interface: 

Currently the RS232 interface is regularly used for 
almost devices on the PC. 

This standard requires a length of 15 m and has a 
low immunity to noise [13]-[14]-[15]. 

5.3.7. The differential interface-RS485: 

These standards provide a universal solution for 
transmission over long distances. It may be taken as 
the solution for our case [1]. 

5.3.8. Reason for the choice of RS485: 

It allows: a transmission distance up to 1200m, 
and a good immunity to electromagnetic 
interference. 

As in our case, the distance is 900m and the site has 
electromagnetic interference, this solution is 
adequate. 

5.3.8. Communication protocol [13]-[15]: 

The link level protocol defines the set of rules on the 
meaning of data in the order in which they are 
transmitted. It prevents access conflicts. 

We can distinguish two communication protocols in 
our case: Sequence of reading the information, and 
the other of writing 

5.3.9. Data processing at Melloussa: 

Once the data is transmitted by the PC, the 
database available at the central unit, will be treated 
later by using the software. Example of software 
Alwin studied in a first part of this work, with its 
system of data acquisition [16]. 

5.4. Industrial solution proposed for remote and 
monitoring the wind farm of Tangier: 

According to industry solutions that we offered for 
the data transmission and control and monitoring of 
wind turbines, we can adopt the following solution 
that encompasses the advantages of the above 
solutions: 

The industrial solution that will apply uses a PLC 



(FLOWTEL) with a management software (XFLOW) to 
secure and maximize the wind farms [ 1 3]-[ 17]: 




Figure 7. Structure of the remote management solution that will be 
applied to the wind farm in Tangier. 



6. Conclusion 

Wind energy is a major contribution towards 
energy independence of Morocco. So, the present 
works consists in feasibility study of harnessing 
wind energy for turbine installation in Tangier’s site 
in northern Morocco. In the first of the study, we 
evaluated: 

• The wind speed potential through the 
treatment of climatologically data of a 
station, using software ALWIN for 
dimensions of wind parks. 

• The prediction of electrical energy 
produced in the site basically on the choice 
of wind energy converter. 

• The determination of the direction of the 
wind energy converter in the site. 

Indeed, we’re studied the integration of a wind 
parks power production in the electric national 
network that have no negative effect on the electric 
grid. However, it empowers the electric national 
network capacity. Eventually, we’re studied and 
suggested some solutions to apply remote 
monitoring of the wind farm to optimize our 
exploitation of the power production. 
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